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Abstract 

Flapping  wing  propulsion  offers  an  attractive  alternative  to  conventional  propul¬ 
sion  methods  for  future  unmanned  aerial  systems.  At  scales  relevant  for  future  micro- 
and  nano-aerial  vehicles,  large  insects  often  utilize  a  single  aerodynamic  surface  under¬ 
going  large  amplitude  oscillatory  motion  to  produce  the  lifting,  propulsive,  and  control 
forces  necessary  for  sustained  flight.  These  dynamic  wing  motions  produce  unsteady 
aerodynamic  phenomena  which  are  responsible  for  the  primary  lifting,  propulsive, 
and  control  forces.  The  use  of  unsteady  aerodynamic  phenomena  differentiate  flap¬ 
ping  wing  fliers  from  conventional  fixed  and  rotary  wing  configurations. 

It  has  been  previously  demonstrated  that  the  unsteady  aerodynamic  phenomena 
are  sensitive  to  variations  in  the  wing  motion  as  well  as  the  dynamic  response  of  flex¬ 
ible  wing  structures.  Yet,  it  is  not  well  understood  how  the  unsteady  aerodynamic 
phenomena,  dynamic  wing  motions,  and  structural  response  interact  to  determine 
overall  aerodynamic  performance.  This  work  examines  the  influence  of  structural  re¬ 
sponse  on  the  local  unsteady  aerodynamic  phenomena  produced  by  the  dynamic  wing 
motions,  and  considers  the  influence  of  local  variations  in  the  unsteady  aerodynamic 
phenomena  on  overall  aerodynamic  performance. 

To  gain  insight  into  the  interplay  between  wing  motion,  structural  response, 
and  unsteady  aerodynamic  phenomena,  Navier-Stokes  based  numerical  simulation 
was  used  to  examine  the  influence  of  prescribed  structural  deformations  on  the  aero¬ 
dynamic  performance  of  a  characteristic  flapping  wing  configuration  based  on  the 
Hawkmoth  Manduca  sexta.  Prescribed  deformations,  including  span-wise  feathering 
and  span-wise  bending,  for  an  elastic  axis  located  at  the  leading  edge  of  the  wing 
were  superimposed  onto  the  rigid  wing  motion  defined  by  a  set  of  experimentally  de¬ 
termined  kinematic  parameters  for  a  Hawkmoth  at  hover.  Both  deformation  modes 
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influence  the  velocity  and  acceleration  profile  of  the  wing  surface,  altering  the  un¬ 
steady  aerodynamic  phenomena  produced  by  the  dynamic  wing  motion. 

The  spanwise  feathering  rotation,  or  torsional  response,  altered  the  motion  of 
the  wing  near  the  wing  root.  This  variation  in  the  acceleration  profile  influenced  the 
non-circulatory  aerodynamic  response  and  the  local  wake  structures  produced  near 
the  wing  root  during  pronation  and  supination.  Increased  lifting  forces  and  enhanced 
aerodynamic  efficiencies  were  observed  for  a  moderate  increase  in  torsional  flexibility. 

Peak  bending  deformations  near  the  wing  tip  also  occured  during  pronation  and 
supination,  altering  the  velocity  and  acceleration  profiles  of  the  wing  as  the  circulatory 
aerodynamic  phenomena  undergo  a  transition  as  the  wing  changes  direction  of  motion. 
Because  of  the  timing  of  the  bending  deformations,  small  tip  deformations  may  have 
a  significant  influence  on  overall  aerodynamic  performance. 
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Influence  of  Structural  Flexibility 
on  Flapping  Wing  Propulsion 

I.  Introduction 

Natural  fliers  operate  effectively  in  flow  regimes  where  the  propulsion  efficiency 
of  conventional  systems  decrease  with  size,  and  demonstrate  an  aerial  agility  that 
has  yet  to  be  achieved  by  engineered  fliers.  One  of  the  features  that  differentiate 
small  flapping  wing  fliers  from  engineered,  fixed  and  rotary  wing,  fliers  is  the  use  of 
unsteady  aerodynamic  phenomena  to  produce  the  aerodynamic  forces  necessary  for 
sustained  flight.  By  understanding  how  natural  fliers  utilize  and  control  the  unsteady 
aerodynamic  phenomena  responsible  for  their  enhanced  aerial  abilities,  engineered 
mini-  and  micro-aerial  systems  may  be  developed  that  take  advantage  of  similar  un¬ 
steady  aerodynamic  mechanisms.  Thus,  bio-inspiration  offers  a  means  to  enhance  the 
performance  of  the  next  generation  of  mini-  and  micro-aerial  systems  over  existing 
fixed  and  rotary  wing  systems. 

1.1  Unsteady  Aerodynamics 

The  unsteady  aerodynamic  phenomena  that  allow  insects  to  operate  efficiently 
at  small  scales  are  produced  by  dynamic  rotations  of  the  insect  wing.  These  wing 
rotations  are  oscillatory  in  nature  and  encompass  a  large  variety  of  motion  profiles 
and  associated  tip  paths  for  different  species  of  insects.  The  tip  path  associated  with 
a  particular  insect  species  depends  on  the  insect  morphology,  the  configuration  of  the 
wing  and  body  structures  and  joints,  and  the  insect  physiology  which  determines  how 
the  wing  is  actuated. 

At  wing  lengths  on  the  order  of  50  millimeters,  the  dominant  unsteady  aerody¬ 
namic  phenomenon  assumed  to  be  responsible  for  aerodynamic  lift  is  a  leading  edge 
vortex  produced  by  a  laminar  flow  separation  near  the  leading  edge  of  the  wing  [1], 
This  vortical  structure  produces  a  region  of  low  pressure  near  the  wing  surface  and 
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influences  the  strength  of  the  bound  circulation  about  the  wing.  This  leading  edge 
vortex  is  similar  to  the  vortical  structure  produced  during  dynamic  stall  observed  for 
conventional  configurations  undergoing  rapid  pitch  or  plunge  maneuvers.  Unlike  the 
vortex  produced  during  a  dynamic  stall,  the  leading  edge  vortex  observed  for  flapping 
wing  fliers  is  stable  for  the  duration  of  the  half-stroke. 

In  addition  to  circulatory  aerodynamic  phenomena,  there  are  aerodynamic  forces 
associated  with  the  acceleration  of  the  volume  of  air  in  the  immediate  vicinity  of  the 
wing  as  the  wing  changes  direction  of  motion  at  the  end  of  each  half-stroke.  This  non- 
circulatory  phenomenon  influences  the  pressure  distribution  along  the  wing  surface 
as  well  as  the  pressure  and  velocity  field  about  the  wing.  The  circulatory  and  non- 
circulatory  aerodynamic  phenomena  each  produce  localized  wake  structures.  These 
local  wake  structures  interact  to  form  the  complex  wake  structures  observed  for  natu¬ 
ral  fliers.  Sustained  flight  is  achieved  when  the  vortical  wake  structures  align  in  space 
and  time  to  form  a  coherent  momentum  jet. 

1.2  Dynamic  Structural  Response  and  Fluid- Structure  Interaction 

Because  the  aerodynamic  phenomena  are  produced  by  the  dynamic  wing  mo¬ 
tion,  they  are  sensitive  to  variations  in  the  motion  profile  [2, 3, 4, 5, 6, 7].  In  addition  to 
variations  in  the  underlying  motion  profile,  any  mechanisms  that  alter  the  wing  mo¬ 
tion,  such  as  structural  response,  may  influence  the  unsteady  aerodynamic  phenomena 
and  the  resulting  aerodynamic  forces  [8,9,10,11].  Most  insects  have  light-weight,  flexi¬ 
ble,  wing  structures  that  exhibit  varying  degrees  of  dynamic  response  during  the  wing 
stroke.  The  amplitude  and  modal  contribution  of  the  dynamic  structural  response 
depend  on  several  factors,  including  the  size  of  the  insect,  the  wing  structure,  material 
properties,  and  the  forces  and  moments  acting  on  the  wing. 

There  are  two  significant  sets  of  forces  and  moments  acting  on  the  wing  struc¬ 
ture,  the  inertial-elastic  forces  due  to  the  dynamic  wing  motion,  and  the  aerodynamic 
forces  produced  in  response  to  the  dynamic  wing  motions.  Both  the  inertial-elastic 
and  aerodynamic  forces  depend  on  the  profile  of  the  motion  at  the  wing  root.  How- 
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ever,  the  aerodynamic  forces  are  also  sensitive  to  phenomena  that  introduce  local 
variations  in  the  wing  motion,  such  as  structural  response.  This  sets  up  a  potential 
feedback  loop  where  aerodynamic  forces  influence  structural  response,  and  structural 
response  influences  the  aerodynamic  forces.  The  influence  of  this  multi-physics  cou¬ 
pling  on  aerodynamic  performance  should  be  understood  in  order  to  optimize  future 
mini-  and  micro-aerial  vehicle  designs  that  utilize  flapping  wing  propulsion. 

1.3  Multi-Physics  Interactions 

Multi-physics  coupling  in  flapping  wing  propulsion  can  be  extended  to  include 
the  kinematic  parameters  that  define  the  dynamic  wing  motion.  If  the  goal  is  to 
optimize  a  design  to  achieve  a  desired  set  of  performance  goals,  to  trim  an  exist¬ 
ing  design  for  a  given  state,  such  as  maximum  efficiency  or  maximum  thrust,  or  to 
achieve  an  agile  maneuver,  the  kinematics,  aerodynamics,  and  structural  dynamics 
should  be  considered  to  be  linked.  Alterations  to  the  actuation  of  the  wing,  or  base 
excitation,  influence  the  unsteady  aerodynamic  phenomena  and  aerodynamic  forces, 
the  inertial-elastic  forces  acting  on  the  wing  structure,  and  the  interactions  between 
the  dynamic  structural  response  and  the  unsteady  aerodynamic  forces.  The  nature  of 
the  fluid-structure-kinematic  coupling  is  characteristic  of  a  specific  insect  species  and 
the  kinematic  parameters  associated  with  the  insect’s  physiology. 

1.4  Multi-Physics  Simulation 

Examining  the  influence  of  multi-physics  coupling  using  numerical  simulation 
depends  on  an  understanding  of  the  influence  of  each  of  the  physical  components 
that  contribute  to  the  coupled  system.  In  the  case  of  the  flapping  wing  problem, 
the  influence  of  kinematics  and  structural  flexibility  on  aerodynamic  performance  are 
active  areas  of  investigation. 

Understanding  these  multi-physics  interactions  requires  a  benchmark  configura¬ 
tion  where  the  characteristics  of  each  of  the  physical  components  are  well  documented. 
For  a  flapping  wing  flier,  this  would  start  with  a  characterization  of  the  wing  structure 
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and  material  properties,  documentation  of  the  kinematic  parameters  that  describe  the 
wing  motion  in  time,  an  understanding  of  the  root  actuation  of  the  wing  structure 
as  a  function  of  the  kinematic  parameters,  and  an  understanding  of  the  aerodynamic 
phenomena  associated  with  the  chosen  wing  kinematics.  Suitable  simulation  tech¬ 
niques  that  have  the  fidelity  to  capture  the  essential  characteristics  of  each  physical 
component  would  then  be  selected  and  coupled  in  order  to  simulate  the  multi-physics 
interactions  and  document  the  influence  of  coupling  on  aerodynamic  performance. 

1.5  Scope  of  Research 

The  objective  of  this  research  is  to  understand  how  structural  response  in  the 
form  of  local  deformations  and  feathering  rotations  influence  the  underlying  unsteady 
aerodynamic  phenomena  that  are  responsible  for  the  aerodynamic  forces  and  influence 
aerodynamic  performance.  This  work  represents  a  benchmark  multi-physics  simula¬ 
tion  of  flapping  wing  propulsion  and  the  first  examination  of  the  influence  of  feathering 
modes  on  the  local  aerodynamic  phenomena  using  three-dimensional  Navier-Stokes 
simulation. 

1.5.1  Thesis  Statement.  The  wings  of  natural  fliers  exhibit  significant  struc¬ 
tural  deformations  in  free  flight.  This  structural  response  influences  the  strength  and 
timing  of  the  local  aerodynamic  mechanisms  produced  by  the  dynamic  wing  motion. 
Variations  in  the  local  aerodynamic  mechanisms  influence  the  overall  aerodynamic 
performance,  as  represented  by  mean  lift,  power,  and  efficiency. 

1.6  Research  Approach 

In  order  to  examine  the  influence  of  structural  response  on  the  unsteady  aero¬ 
dynamic  phenomena  and  the  associated  influence  on  the  aerodynamic  forces  and  mo¬ 
ments,  numerical  simulation  of  the  Navier-Stokes  equations  was  used  to  predict  the 
aerodynamic  performance  of  a  test  configuration  based  on  the  Hawkmoth  Manduca 
sexta  with  prescribed  structural  deformations.  The  available  body  of  literature  ex- 
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amining  the  Hawkmoth  include  characterizations  of  the  unsteady  aerodynamic  mech¬ 
anisms,  kinematic  parameters  that  approximate  the  rigid  wing  motion  in  the  flight 
regime  from  hover  to  forward  flight,  and  static  characterizations  of  the  Hawkmoth 
wing  structure.  This  body  of  literature,  while  incomplete  from  the  standpoint  of  a 
coupled  multi-physics  simulation,  allows  for  the  influence  of  the  dominant  structural 
mode,  feathering,  or  torsion  about  the  wing  length,  on  the  unsteady  aerodynamic 
mechanisms  to  be  examined  using  a  decoupled  simulation  approach. 

As  an  initial  investigation,  several  preliminary  tasks  had  to  be  accomplished  in 
order  to  examine  the  influence  of  structural  response  on  the  unsteady  aerodynamic 
phenomena.  The  preliminary  tasks  fell  into  two  main  categories,  identification  and 
development  of  a  suitable  test  configuration  based  on  available  literature,  and  the 
identification,  acquisition,  modification,  and  demonstration  of  a  suitable  flow  solver. 

To  examine  the  influence  of  structural  response  on  the  aerodynamic  mechanisms, 
it  was  considered  prudent  by  the  author  to  understand  the  influence  of  kinematic  ap¬ 
proximation  on  the  aerodynamic  forces.  A  complete  investigation  considering  the 
influence  of  individual  kinematic  parameters  was  too  intensive  for  a  study  also  con¬ 
sidering  the  influence  of  structural  flexibility.  In  order  to  gain  insight  into  the  influence 
of  the  dynamic  wing  motions  on  the  unsteady  aerodynamic  phenomena,  the  kinematic 
approximation  used  for  this  study  was  systematically  varied  for  a  selected  set  of  sim¬ 
ulations.  Simulations  were  computed  for  wing  motion  based  on  2-,  3-,  and  4-term 
Fourier  approximations  of  experimentally  derived  kinematic  parameters.  These  rigid- 
wing  simulations  also  served  as  a  baseline  when  examining  the  influence  of  structural 
response  on  the  unsteady  aerodynamic  phenomena  and  aerodynamic  performance. 

While  a  complete  physical  characterization  of  the  wing  structure  was  not  avail¬ 
able  for  this  research,  the  published  study  documenting  the  dynamic  wing  motions 
at  hover  recorded  feathering  orientation  along  the  wing  length  as  a  function  of  time. 
This  measurement  allowed  a  spanwise  feathering  distribution  along  the  wing  length 
to  be  reconstructed.  This  reconstruction  of  the  feathering  parameter  was  used  to  pre- 
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scribe  the  feathering  rotation  of  the  wing  surface  for  the  numerical  simulations.  The 
influence  of  bending  was  examined  using  an  inertially  motivated  bending  distribution 
and  time  history,  and  is  presented  as  a  proof-of-concept  demonstration  of  the  solution 
method  and  as  an  extension  of  the  fluid-structure  coupling  mechanisms  studied  with 
the  prescribed  torsional  deformations. 

1 . 7  Document  Organization 

The  remainder  of  this  document  is  organized  as  follows: 

Chapter  II:  Covers  the  fundamentals  of  insect  flight  and  the  relevant  aerody¬ 
namic  phenomena.  This  chapter  also  provides  an  overview  of  the  aerodynamic  sim¬ 
ulation  method  and  the  capabilities  of  OVERFLOW,  the  Navier-Stokes  flow  solver 
selected  for  this  research,  as  well  as  the  modifications  that  form  the  basis  of  the 
updated  flow  solver,  OVERFLOW  2.1-Elastic. 

Chapter  III:  Provides  a  discussion  of  the  verification  and  validation  of  the  up¬ 
dated  OVERFLOW  2.1-Elastic  flow  solver,  and  demonstrates  the  ability  of  OVER¬ 
FLOW  2.1-Elastic  to  capture  the  characteristic  aerodynamics  exhibited  by  a  bench¬ 
mark  configuration  modeled  on  the  Hawkmoth. 

Chapter  IV:  Details  the  method  in  which  the  complex  behavior  of  the  Hawk- 
moth  was  reduced  to  a  characteristic  configuration  that  exhibits  the  unsteady  aero¬ 
dynamic  phenomena  observed  in  previous  studies  while  providing  a  simplified  config¬ 
uration  for  future  numerical  and  experimental  studies. 

Chapter  V:  Examines  the  influence  of  the  kinematic  approximation  on  the 
unsteady  aerodynamic  mechanisms  and  aerodynamic  performance. 

Chapter  VI:  Examines  the  influence  of  span-wise  feathering  variation  on  the 
unsteady  aerodynamic  phenomena. 

Chapter  VII:  Presents  a  proof-of-concept  examination  of  the  influence  of  bend¬ 
ing  and  the  influence  on  the  unsteady  aerodynamic  phenomena. 
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Chapter  VIII:  Summarizes  the  influence  of  the  kinematic  profile,  feathering 
variation,  and  bending  on  the  unsteady  aerodynamic  mechanisms  and  presents  rec¬ 
ommendations  for  future  research  efforts. 
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II.  Numerical  Simulation  of  a  Flexible  Flapping  Wing  Flier 

Natural  fliers  range  in  size  from  the  millimeter  scale  to  the  meter  scale  and  operate  at 
Reynolds  numbers  from  100  to  100,000.  Throughout  this  range  of  sizes  and  operating 
conditions,  natural  fliers  utilize  a  variety  of  wing  motions  to  generate  the  aerodynamic 
forces  necessary  for  sustained  flight.  The  selection  of  a  single  representative  species  for 
a  detailed  examination  represents  a  balancing  act  between  desired  characteristics  and 
the  availability  of  previous  research  upon  which  a  numerical  model  may  be  developed. 

2. 1  The  Hawkmoth 

The  Hawkmoth  Manduca  sexta  represents  one  such  species  that  exhibits  capa¬ 
bilities  desired  for  the  next  generation  of  mini-  and  micro-aerial  vehicles  at  a  scale 
relevant  for  future  micro-aerial  vehicle  designs.  A  typical  adult  Hawkmoth  has  a 
wing  span  around  110  millimeters  and  total  mass  of  1.8  grams  [12,13],  which  falls 
within  the  upper  bound  of  the  DARPA  definition  of  a  micro-aerial  vehicle  [14].  The 
aerodynamic  phenomena  utilized  by  the  Hawkmoth  to  achieve  sustained  hover,  agile 
maneuvers,  and  forward  flight  have  been  previously  studied  by  Ellington,  Liu,  and 
others  [15,16,17,6,18,19,20,21,22,23,24,25],  who  have  documented  the  dynamic  wing 
motions,  unsteady  aerodynamic  phenomena,  and  morphological  characteristics. 

2.1.1  Aerodynamic  Characteristics.  In  this  operating  regime,  one  of  the 
dominant  unsteady  aerodynamic  mechanisms  is  the  leading  edge  vortex  produced 
during  the  translational  phase  of  the  wing  stroke  [15].  The  leading  edge  vortex  is 
produced  by  a  laminar  flow  separation  just  aft  of  the  leading  edge  that  forms  a  strong 
vortical  structure  as  shown  in  Figure  1.  This  vortical  structure  produces  a  region  of 
low  static  pressure  that  influences  the  pressure  distribution  along  the  wing,  shown 
in  Figure  2,  due  to  the  proximity  of  the  vortical  structure  to  the  wing  surface.  As 
previously  discussed  in  Section  1.1,  the  leading  edge  vortex  is  similar  to  dynamic 
stall  observed  for  conventional  configurations  at  high  angles  of  attack.  However, 
the  leading  edge  vortex  exhibited  by  the  Hawkmoth  is  stabilized  for  the  duration  of 


Figure  1:  Diagram  of  a  typical  leading  edge  vortex.  Vortical  structures  are  shown 

using  blue  dashed  lines  and  bounded  by  green  solid  lines.  Streamlines  are  shown  using 
solid  blue  lines,  and  the  wing  surface  is  shown  in  gray.  Figure  adapted  from  van  den 
Berg  and  Ellington  [16]. 

each  half-stroke  by  a  span-wise  flow  towards  the  wing  tip  at  the  core  of  the  vortical 
structure. 

The  finite  length  Hawkmoth  wing  exhibits  significant  three-dimensional  tip  ef¬ 
fects  which  produce  a  strong  vortical  structure  at  the  wing  tip.  The  three-dimensional 
tip  effects  are  accentuated  by  the  spanwise  pressure  distribution,  which  features  a  re¬ 
gion  of  low  pressure  near  the  wing  tip,  as  shown  in  Figure  2.  The  spanwise  flow  and 
the  leading  edge  vortex  coalesce  with  the  tip  vortex  to  form  a  continuous  vortical 
structure  that  extends  along  the  length  of  the  wing  from  near  the  wing  root  until  it 
is  shed  into  the  wake  along  the  path  of  the  wing  tip,  also  shown  in  Figure  1. 

The  end  of  each  half-stroke  is  characterized  by  a  deceleration  of  the  wing  sweep 
until  the  translational  velocity  component  reaches  zero  and  the  direction  of  motion 
changes.  This  deceleration  is  accompanied  by  a  feathering  rotation  along  a  span-wise 
axis  that  re-orients  the  leading  edge  of  the  wing  for  the  following  half-stroke.  As 
the  wing  changes  direction  of  motion,  the  strength  of  the  leading  edge  vortex  that 
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Figure  2:  Pressure  distribution  along  a  Manduca  sexta  planform  computed  using 

rigid  wing  motions. 
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Figure  3:  Evolution  of  the  leading  edge  vortex,  starting  vortex,  and  stopping  vortex 
during  the  stroke  cycle.  Vortical  structures  produced  during  the  upstroke  are  shown 
as  solid  lines  and  vortical  structures  produced  during  the  downstroke  are  shown  as 
dashed  lines.  The  velocity  of  the  wing  at  the  leading  and  trailing  edges  is  shown  in 
red.  Figure  adapted  from  van  den  Berg  and  Ellington  [16]. 


developed  during  the  preceding  half-stroke  decreases  until  it  is  released  into  the  wake. 
This  decrease  in  the  strength  of  the  leading  edge  vortex  and  bound  circulation  about 
the  wing  is  balanced  by  a  release  of  vorticity  into  the  wake  at  the  trailing  edge.  The 
vorticity  shed  from  the  trailing  edge  of  the  wing  forms  a  stopping  vortex,  shown  in 
Figure  3.  As  the  direction  of  motion  changes  and  the  wing  accelerates  into  the  next 
half-stroke,  a  new  leading  edge  vortex  is  formed  along  with  an  accompanying  starting 
vortex  that  balances  the  change  in  bound  circulation  about  the  wing.  This  cycle  is 
illustrated  in  Figure  3. 

The  vortex  wake  produced  by  the  leading  edge  vortex,  tip  vortex,  and  the 
stopping  and  starting  vortices  produced  by  each  wing  form  a  vortex  ring  [16,  18, 
21].  When  the  phase  and  orientation  of  the  vortical  structures  align,  the  vortex  ring 
produces  a  momentum  jet  that  provides  the  thrust  necessary  for  sustained  hovering 
and  forward  flight,  as  shown  by  the  streamlines  outside  the  vortex  system  in  Figure 
1. 
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In  addition  to  the  circulatory  aerodynamic  mechanisms  that  produce  the  lead¬ 
ing  edge  vortex  and  influence  the  bound  circulation  about  the  wing,  the  dynamic 
wing  motions  produce  additional  aerodynamic  forces  associated  with  the  acceleration 
and  deceleration  of  the  volume  of  air  immediately  surrounding  the  wing.  These  non- 
circulatory  forces  influence  the  pressure  distribution  along  the  wing  length  as  well  as 
the  pressure  distribution  in  the  fluid  immediately  surrounding  the  wing.  These  pres¬ 
sure  variations  influence  the  local  velocity  distribution,  and  potentially  alter  the  for¬ 
mation  and  convection  of  the  vortical  wake  structures  associated  with  the  circulatory 
aerodynamic  phenomena.  Wing-wake  interaction  is  an  example  of  a  non-circulatory 
aerodynamic  phenomenon  that  occurs  at  the  end  of  each  half-stroke  as  the  wing  “en¬ 
counters”  the  wake  that  is  entrained  behind  the  wing  during  the  half-stroke. 

2.1.2  Structural  Characteristics.  Hawkmoth  wings  are  lightweight,  compli¬ 
ant  structures  that  exhibit  significant  structural  response  during  the  wing  stroke.  A 
typical  wing  response  at  hover  is  captured  in  the  sequence  of  stills  shown  in  Figures  4 
and  5.  The  sequence  of  still  frames  highlight  the  large  variation  in  feathering  orienta¬ 
tion  of  the  wing  as  a  function  of  span-wise  position,  while  the  shape  and  magnitude  of 
any  deformations  associated  with  bending  is  difficult  to  discern  through  observation 
alone.  Figure  6  provides  an  example  of  how  point  of  view  can  influence  qualitative 
estimates  of  bending. 

Previous  experimental  investigations  indicate  that  the  Hawkmoth  wing  is  a 
non-homogeneous  structure  with  pronounced  anisotropic  properties  [26],  with  sev¬ 
eral  orders  of  magnitude  difference  between  the  static  bending  stiffnesses  measured  in 
the  chord-wise  and  span-wise  directions.  Lacking  a  dynamic  characterization  of  the 
Hawkmoth  wing  structure,  trends  must  be  drawn  from  structural  characterizations 
that  have  been  performed  on  other  insects,  such  as  the  Fruit  Fly,  Dragon  Fly,  and 
the  Desert  Locust  [27,28,29,30,31,32,33,34,35],  and  bulk  structural  properties  for 
the  Hawkmoth  wing  estimated  using  static  bending  and  torsion  tests  [36,37,38]. 
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Figure  4:  A  typical  upstroke  for  a  Hawkmoth  at  hover.  The  point  of  view  highlights 
the  span-wise  feathering  variation  at  t/r  =  1/6,  and  the  wing  planform  at  t/r  =  3/6. 
Video  stills  are  reproduced  with  permission  [13]. 
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Figure  5:  A  typical  downstroke  for  a  Hawkmoth  at  hover.  The  point  of  view  hinders 
qualitative  observation  of  the  feathering  variation  during  the  half-stroke.  Video  stills 
are  reproduced  with  permission  [13]. 


14 


Frame  6  of  42 

•41*; 

•  /  h.n- 

/ 

Mk  .  \ 

/  ■ 

f  \ 

Curvature  due  to  bending  or  planform? 

Frame  18  of  42 

/  /  g 

/  -lib  t 

Leading  edge  curvature  due  to  planform 

/ 

Wing  1  parting  Frlgp 

Chordwise  Section  - 

Figure  6:  Estimating  the  relative  magnitude  of  bending  and  torsion  modes  dur¬ 

ing  the  wing  stroke  is  difficult  to  accomplish  through  qualitative  observation  alone. 
Choosing  a  single  point  of  view  to  observe  the  dynamic  structural  response  is  compli¬ 
cated  by  the  large  amplitude  wing  rotations,  and  determining  the  relative  contribu¬ 
tions  of  the  bending  and  torsion  modes  is  complicated  by  the  curvature  of  the  leading 
edge  of  the  Hawkmoth  wing  when  using  multiple  points  of  view. 


15 


2.2  Analytic  Aerodynamic  Models 

A  variety  of  analytic  models  have  been  proposed  for  examining  insect  flight.  The 
use  of  these  methods  has  evolved  along  with  our  understanding  of  the  aerodynamic 
phenomena  associated  with  flapping  wing  flight.  Early  aerodynamic  analytic  models 
were  based  on  quasi-steady  aerodynamic  assumptions.  The  blade  element  model  is 
one  such  technique  used  to  examine  flapping  wing  fliers  that  idealized  the  wings  as 
a  set  of  “blade  elements”,  or  chordwise  sections  [39].  Aerodynamic  forces  for  each 
chordwise  section  were  based  on  a  quasi-steady  analysis  of  each  blade  element  based 
on  wing  orientation  and  velocity.  The  wing  velocity  represents  a  combination  of  the 
flapping,  forward,  and  induced  velocity  components.  The  induced  velocity  component 
were  often  based  on  a  momentum  jet  analysis  of  the  wake. 

Models  based  on  a  momentum  jet  analysis  allowed  for  mean  aerodynamic  per¬ 
formance  and  power  requirements  to  be  examined.  The  momentum  jet  approach 
idealizes  the  wing  as  an  actuator  disc  which  imparts  momentum  from  the  wing  to 
the  fluid.  Aerodynamic  forces  and  power  may  be  predicted  from  an  analysis  of  the 
farfield  flow  velocities.  Because  the  wing  is  idealized  as  an  actuator  disc  and  depends 
on  farfield  wake  characteristics,  the  influence  of  kinematic  and  geometric  variations 
on  the  near-body  wake  are  not  included  in  the  analysis  [39]. 

Hybrid  methods,  which  combine  the  blade  element  method  with  a  modified  mo¬ 
mentum  theory  approach  were  proposed  that  allowed  the  influence  of  morphological, 
kinematic,  and  geometric  variations  on  the  mean  aerodynamic  force  and  power.  The 
hybrid  momentum  jet  was  a  modification  of  the  blade  element  model  which  utilized 
a  pulsed  actuator  disc  to  predict  the  induced  velocity  component  [40,41], 

Other  efforts  were  made  to  examine  the  time-dependent  aerodynamic  forces, 
such  as  the  unsteady  strip  theory  based  on  the  Theodorsen  function  for  harmonic 
pitching  and  plunging  motions  described  by  DeLaurier  [42],  and  three-dimensional 
influences,  such  as  the  vortex  lattice  method  described  by  Sunada  et  al.  [43,44]  to  ex¬ 
amine  three-dimensional  oscillating  flat  plates.  The  strip  method  described  by  DeLau- 
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rier  overcame  the  quasi-steady  aerodynamic  assumption  by  utilizing  the  Theodorsen 
function  for  harmonic  pitching  and  plunging  airfoils.  However,  the  Theodorsen  ap¬ 
proximation  is  only  valid  for  small  oscillation  amplitudes. 

2.3  Physics  Based  Numerical  Simulation 

Physics  based  simulation  is  an  alternative  to  analytical  methods  for  examining 
the  aerodynamic  performance  of  flapping  wing  fliers.  Smith  [45]  described  the  use  of 
an  unsteady  panel  method  to  examine  the  performance  of  a  Hawkmoth  during  flight. 
The  unsteady  panel  method  utilized  doublet  panels  to  model  the  wing  surface  and  the 
evolution  of  the  unsteady  wake  in  time.  The  unsteady  panel  method  overcomes  many 
of  the  limitation  of  the  analytic  models  by  treating  the  aerodynamic  phenomenon 
as  unsteady  and  allowing  for  a  free-wake  analysis.  However,  vortex  panel  methods 
are  based  on  inviscid  aerodynamic  approximations  and  are  unable  to  capture  viscous 
effects,  such  as  the  laminar  flow  separation  that  produces  the  leading  edge  vortex 
characteristic  of  many  natural  fliers  at  low  Reynolds  numbers. 

Time-accurate  Navier-Stokes  based  simulation  of  flapping  wing  fliers  have  been 
demonstrated  using  a  variety  of  approaches  over  the  course  of  the  last  decade.  Finite 
volume  solutions  using  body  conformal  grids  have  been  computed  by  Liu  [20],  Liu  et 
al.  [19],  Aono  [24],  and  Aono  et  al.  [25],  Sun  and  Tang  [46,47]  and  Sun  et  al.  [48], 
among  others.  Other  formulations  of  the  Navier-Stokes  equations  have  been  utilized 
for  simulating  natural  fliers,  such  as  the  finite  element  formulation  demonstrated  by 
Ramamurti  and  Sandberg  [49,50]  and  Hamamoto  [51],  and  finite  difference  solutions 
using  non-body  conforming  solution  domains  and  immersed  boundary  conditions  by 
Mittal  [52,53]  and  Mittal  et  al.  [54], 

2-4  Overview  of  the  Aerodynamic  Simulation  Method 

For  the  purpose  of  numerical  simulation,  the  large  amplitude  wing  motions  may 
be  represented,  or  resolved,  into  two  components;  a  rigid  body  rotation  with  respect 
to  the  wing  root,  and  local  deformations  and  rotations  with  respect  to  the  undeformed 
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wing.  The  large  amplitude  rigid  body  motions  for  three-dimensional  configurations 
have  been  accommodated  using  a  variety  of  techniques,  including  mesh  regeneration, 
mesh  deformation,  overset  methods,  and  immersed  boundary  techniques.  Each  tech¬ 
nique  has  advantages  and  disadvantages  in  terms  of  gridding  the  complex  geometric 
configurations,  solution  efficiency,  and  the  extensibility  of  the  method  for  including 
fluid-structure  interaction  and  modeling  multiple  moving  components. 

After  reviewing  the  solution  methods  and  solvers  utilized  in  previous  investiga¬ 
tions  and  the  available  commercial  flow  solvers,  the  OVERFLOW  solver  was  selected 
for  use  in  the  current  research  effort.  The  OVERFLOW  solver  was  selected  because 
it  offers  time-accurate  low  Mach  number  preconditioning  [55],  higher-order  solution 
capabilities  [56],  dynamic  rigid  body  solution  capabilities  [57],  and  aeroelastic  simula¬ 
tion  capabilities  [58].  In  addition  to  the  solver  capabilities,  OVERFLOW  is  a  NASA 
research  code.  As  such,  the  use  of  the  software  is  not  license-limited  in  the  AFIT  HPC 
or  the  AFRL  MSRC  environments,  and  the  source  code  is  available  for  examination 
and  modification. 

2.5  Overview  of  the  OVERFLOW  Flow  Solver 

The  OVERFLOW  code  is  a  compressible  Navier-Stokes  flow  solver  for  struc¬ 
tured  overset  solution  domains  based  on  a  finite  difference  Navier-Stokes  formulation. 
The  OVERFLOW  solver  has  undergone  multiple  revisions  and  branching  since  the 
original  code  development,  as  illustrated  in  Figure  7.  The  versions  of  the  code  of  par¬ 
ticular  interest  for  the  current  study  are  OVERFLOW-D,  which  introduced  dynamic 
simulation  capabilities  for  rigid  moving  bodies,  OVERFLOW-D2  which  introduced 
fluid-structure  interaction  capabilities,  and  the  current  development  branch,  OVER¬ 
FLOW  2.1,  which  introduced  higher-order  solution  capabilities. 

The  current  development  version  of  OVERFLOW,  version  2.1,  is  capable  of 
fifth  order  spatial  accuracy  using  the  Weighted  Essentially  Non-Oscillatory  (WENO) 
scheme  [56]  and  second-order  temporal  accuracy  using  Dual  Time  Stepping  [59]. 
OVERFLOW  2.1  is  capable  of  modeling  flow-fields  characterized  by  mixed  Mach 
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numbers,  such  as  those  experienced  by  a  rotorcraft  in  hover,  through  the  use  of  low 
Mach  number  preconditioning.  In  addition  to  higher-order  solutions  at  low  or  mixed 
Mach  numbers,  OVERFLOW  2.1  is  capable  of  modeling  dynamic  rigid-body  motion 
for  multi-body  problems  with  relative  motion,  including  six-degree-of-freedom  sim¬ 
ulations  for  motion  due  to  aerodynamic  forcing.  A  brief  description  of  the  relevant 
OVERFLOW  features  for  modeling  flapping  wing  flight  will  be  provided  here,  while 
a  full  examination  of  the  OVERFLOW  solver  may  be  found  in  the  OVERFLOW  2.1 
documentation  [60,61]  and  related  publications  [62,63,64,65,66,55,67,68]. 

2.5.1  Chimera  Solution  Approach.  The  overset,  or  Chimera  solution  meth¬ 
ods  enable  structured  flow  solvers  to  compute  flow  solutions  for  complex  domains. 
Overset  techniques  simplify  the  generation  of  solution  domains  about  complex  ge¬ 
ometries  and  configurations  by  allowing  separate  computational  domains  to  overlap 
in  physical  space.  By  allowing  computational  domains  to  overlap  in  physical  space, 
complex  geometric  features,  or  combinations  of  features  that  are  difficult  to  model 
using  a  single  structured  solution  domain  may  be  modeled  using  multiple  overlapping 
computational  domains  [69]. 

Because  solution  domains  are  allowed  to  overlap  in  the  physical  space,  an  extra 
boundary  condition  is  required  to  accommodate  the  overlapping  domains.  The  addi¬ 
tional  boundary  condition  is  a  variation  on  a  block-to-block  boundary  condition  for 
multi-block  solvers,  but  instead  of  passing  information  across  an  edge  or  face,  the  over¬ 
set  boundary  condition  interpolates  the  flow  solution  from  points  in  the  overlapping 
solution  regions  [70]. 

Much  of  the  additional  computational  overhead  and  complexity  traditionally 
associated  with  the  Chimera  solution  approach  is  related  to  establishing  interpolation 
stencils  between  solution  domains  for  the  overlapping  boundary  condition,  as  well  as 
the  identification  and  blanking  of  cells,  or  points,  that  fall  within  a  specified  distance 
of  solid  bodies  or  domain  boundaries  in  the  overlapping  regions.  The  Chimera  solution 
approach  requires  that  additional  care  be  taken  when  integrating  aerodynamic  forces 
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Figure  7:  Evolution  the  OVERFLOW  2.1-Elastic  flow  solver  developed  for  this  study.  OVERFLOW  development  history 
based  on  Slotnick  [67]. 


and  moments  on  surfaces  with  overlapping  regions  such  that  the  overlapping  regions 
are  not  integrated  multiple  times  [71]. 

2.5.2  Dynamic  Multi-body  Simulation.  The  Chimera  solution  approach 
is  well  suited  for  computing  dynamic  or  moving  body  flow  simulations.  Using  the 
Chimera  approach,  computational  domains  may  be  grouped  so  as  to  represent  a  sin¬ 
gle  water  tight  body.  These  grouped  solution  domains  may  be  translated  and  rotated 
relative  to  other  moving  “bodies”  or  fixed  solution  domains.  Because  the  solution 
domains  are  dynamic,  the  interpolation  stencils  and  cell  blanking  required  for  the 
overlapping  boundary  condition  must  be  updated  at  each  time  step  for  the  new  solu¬ 
tion  configuration.  This  introduces  additional  computational  overhead  at  each  time 
step.  A  constant  cell  quality  may  be  maintained,  however,  and  the  overset  method 
does  not  require  the  use  of  mesh  deformation  or  mesh  regeneration  to  represent  rigid 
body  motions  [72], 

OVERFLOW  utilizes  the  Domain  Connectivity  Function  utility  (DCF)  [60]  in¬ 
corporated  from  the  OVERFLOW-D  family  of  solvers  to  update  hole-cutting  and 
fringe  point  interpolations  between  near-body  and  off-body  grids  [73].  Arbitrary  rigid 
body  motion  and  domain  grouping  may  be  prescribed  using  the  Geometry  Manipula¬ 
tion  Protocol  (GMP)  [74]  or  through  a  user  defined  subroutine  [60].  OVERFLOW  2.1 
also  includes  a  six-degree  of  freedom  (6DOF)  model  for  computing  rigid  body  motion 
due  to  aerodynamic  forces. 

2.5.3  Low  Mach  Number  Preconditioning.  The  use  of  a  compressible  formu¬ 
lation  of  the  Navier-Stokes  equations,  Equation  (1),  to  compute  the  flow  solution  for 
configurations  with  low  characteristic  free-stream  velocities  generally  results  in  poor 
solution  accuracy  and  solution  convergence  [55].  The  poor  performance  of  the  com¬ 
pressible  formulation  of  the  Navier-Stokes  equations  is  caused  by  a  large  separation 
of  the  eigenvalues,  Equation  (2),  which  depend  on  the  free-stream,  or  convective  ve¬ 
locity,  and  the  acoustic  velocity,  or  speed-of-sound,  which  impose  significant  stability 
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limitations  for  time  marching  solution  methods. 


dq  dE  dF  dG 
dt  drj  d( 
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The  efficient  solution  of  low  Mach  number  configurations  may  be  achieved  by  using 
an  incompressible  formulation  of  the  Navier-Stokes  equations,  such  as  the  Marker- 
and-Cell  (MAC)  method  of  Harlow  and  Welch  [75]  or  the  Semi-Implicit  Method  for 
Pressure  Linked  Equations  (SIMPLE)  method  of  Patankar  [76].  Alternatively,  the 
compressible  formulation  of  the  Navier-Stokes  equations  may  be  modified  with  the 
addition  of  a  pseudo-compressibility  term  that  re-couples  the  conservation  of  mass 
and  conservation  of  momentum  equations  for  incompressible  flows. 


In  order  to  solve  low  Mach  number  and  mixed  Mach  number  flows,  OVERFLOW 
2.1  incorporates  the  Smith- Weiss  preconditioning  scheme  [55,61].  The  Smith- Weiss 
preconditioning  scheme  introduces  a  preconditioning  matrix,  Tp,  acting  on  the  pre¬ 
conditioned  conservative  variables,  qp.  This  formulation  introduces  a  pseudo-time 
variable,  r,  into  the  preconditioned  Navier-Stokes  equation,  Equation  (3). 
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The  preconditioned  Navier-Stokes  equation  produces  a  set  of  eigenvalues,  Equation 
(4),  that  are  balanced  at  low  Mach  numbers. 
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0.5U  (0  +  1)  —  0.5  (P  —  l)2  +  /3c2 
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(4) 


This  formulation  depends  on  a  “preconditioned”  Mach  number,  /?,  which  is  defined 
in  Equation  (5). 

P  =  max  [min  (M2,  l)  ,  Pmin]  (5) 

In  OVERFLOW  2.1,  the  default  value  of  Pmin  is  3 M\EF.  When  used  within  a  dual 
time  stepping  solution  scheme,  the  implemented  Smith- Weiss  preconditioning  pro¬ 
vides  a  time-accurate  solution  when  the  pseudo-time  variables  reach  convergence  at 
each  time  step  [60,61]. 


2.5.4  Aeroelastic  Framework.  The  OVERFLOW-D2  solver  has  been  used  to 
examine  rotorcraft  aeromechanics  using  both  loose  and  tight  coupled  fluid-structure 
interaction  (FSI)  solution  methods  [57,  58,  77].  Loose  aeroelastic  coupling  is  im¬ 
plemented  on  a  per-rotor  revolution  basis,  while  tight  aeroelastic  coupling  may  be 
achieved  at  the  time  iteration  level.  Both  coupling  methods  utilize  the  same  FSI 
framework  which  depends  on  external  programs  to  provide  structural  deformations. 
Fortunately,  the  per-revolution  loose  coupling  method  has  a  corollary  in  the  flap¬ 
ping  wing  problem.  For  a  flapping  wing  flier,  the  wing  kinematics  for  steady  flight 
modes,  both  hover  and  forward  flight,  are  time-periodic  and  may  be  mapped  to  a  rotor 
revolution.  Thus  the  loose  coupling  framework  may  be  used  to  prescribe  structural 
deformations  on  a  per-flapping  cycle  basis. 
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2.5.5  Mesh  Deformation.  The  0VERFL0W-D2  mesh  deformation  capa¬ 
bilities  were  developed  for  slender,  elastic  rotor  blades  exhibiting  “beam”  type  defor¬ 
mations.  Surface  deformations  and  rotations  are  prescribed  along  an  elastic  axis  as 
a  function  of  span,  r  / 11.  and  azimuth,  0  <  T  <  27T.  Because  of  the  focus  on  “beam” 
type  deformations,  the  current  mesh  deformation  implementation  is  unable  to  account 
for  chord-wise  deformations. 

Surface  deformations  are  interpolated  for  the  rotor  azimuth  at  each  time  step 
using  a  cubic  spline  which  maintains  C 1  continuity  of  the  motion  derivative  in  time. 
The  computational  surface  is  deformed  using  a  transformation  matrix  incorporating 
both  transformations  and  rotations.  The  computational  volume  is  then  deformed 
using  the  transformation  matrix  for  each  corresponding  span-wise  location.  This 
algebraic  method  preserves  grid  quality  and  initial  grid  orthogonality  at  the  boundary 
surface  [58]. 

2.6  OVERFLOW  2.1 -Elastic 

In  order  to  utilize  the  fluid-structure  interaction  capabilities  of  OVERFLOW- 
D2  and  the  higher-order  solutions  of  OVERFLOW-2. 1,  the  fluid-structure  interaction 
framework  and  mesh  deformation  routines  from  OVERFLOW-D2,  specifically  OVER¬ 
FLOW  2.0aa  elastic  1.8,  were  incorporated  into  OVERFLOW-2. lh.  This  merged 
solver,  developed  by  the  author  specifically  for  this  research,  is  referred  to  as  OVER¬ 
FLOW  2. 1-Elastic  for  the  remainder  of  this  document. 

2. 7  Summary 

The  Hawkmoth  Manduca  sexta  was  selected  as  the  focus  of  this  investigation 
based  on  the  aerial  capabilities  of  the  Hawkmoth  and  the  availability  of  an  existing 
body  of  literature  documenting  the  unsteady  aerodynamic  characteristics  and  the  dy¬ 
namic  wing  motions  across  a  range  of  flight  conditions.  The  aerodynamic  phenomena 
utilized  by  the  Hawkmoth  to  produce  the  aerodynamic  forces  are  unsteady  in  nature 
and  exhibit  significant  non-linearities  in  response  to  the  dynamic  wing  motion. 
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In  order  to  examine  the  influence  of  structural  response  on  the  unsteady  aerody¬ 
namic  phenomena  and  aerodynamic  performance,  the  aerodynamic  response  will  be 
simulated  using  Navier-Stokes  based  simulation  using  the  OVERFLOW  2.1-Elastic 
flow  solver.  OVERFLOW  2. 1-Elastic  was  developed  specifically  for  this  study,  and 
represents  an  enhancement  of  OVERFLOW  2.1,  which  incorporates  higher  order  solu¬ 
tion  capabilities,  with  the  fluid-structure  interaction  capabilities  of  OVERFLOW-D2. 
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III.  Verification  and  Validation  of  OVERFLOW  2.1-Elastic 


Verification  and  validation  of  OVERFLOW  2.1-Elastic  focused  on  two  main  tasks, 
demonstrating  that  the  modified  solver  retained  the  baseline  solution  capabilities  of 
OVERFLOW  2.1,  and  demonstrating  the  ability  of  OVERFLOW  2.1-Elastic  to  ade¬ 
quately  capture  the  unsteady  aerodynamic  phenomena  exhibited  by  the  Hawkmoth  in 
flight,  as  discussed  in  the  previous  chapter.  The  ability  of  OVERFLOW  2.1-Elastic 
to  capture  unsteady  low-Mach  number  flows  was  demonstrated  by  computing  the 
flow  held  about  a  two-dimensional  circular  cylinder  in  a  cross  how  across  a  range 
of  Reynolds  numbers.  The  circular  cylinder  solutions  demonstrate  the  steady  state, 
time-accurate,  and  dynamic  solution  capabilities  of  OVERFLOW  2.1-Elastic. 

An  additional  demonstration  of  the  dynamic  solution  capabilities  at  low  Mach 
numbers  was  based  on  a  two-dimensional  heaving  airfoil  solution  for  which  experi¬ 
mental  how  measurements  were  available.  The  ability  of  OVERFLOW  2.1-Elastic  to 
capture  the  unsteady  aerodynamic  phenomena  exhibited  by  the  Hawkmoth  in  flight 
was  demonstrated  for  a  baseline  three-dimensional  test  case  that  produced  sufficient 
vertical  force  to  sustain  the  weight  of  a  typical  Hawkmoth  at  hover. 

3.1  Cylinder  in  Cross-flow 

The  circular  cylinder  in  cross-how  is  a  useful  test  case  to  examine  both  the  steady 
and  unsteady  solution  capabilities  of  OVERFLOW  2.1-Elastic.  The  circular  cylinder 
in  a  cross-how  exhibits  distinct  how  characteristics  that  depend  on  Reynolds  number. 
Of  interest  for  this  study  are  the  steady  how  regime  between  Reynolds  number  5  and 
45  in  which  the  wake  forms  a  laminar  separation  bubble  on  the  downstream  side  of  the 
cylinder,  and  the  unsteady  how  regime  between  Reynolds  number  50  and  185  in  which 
the  laminar  wake  forms  a  two-dimensional  von  Karman  vortex  street  downstream  of 
a  three-dimensional  cylinder  [78,79]. 

3.1.1  Laminar  Separation  Bubble.  At  Reynolds  numbers  between  5  and 
45  a  laminar  separation  bubble,  or  recirculation,  forms  on  the  downstream  side  of  a 


26 


circular  cylinder  in  cross-flow.  The  length  of  the  recirculation  region,  L,  relative  to  the 
diameter,  D ,  of  the  cylinder,  L/D,  depends  on  the  diameter  based  Reynolds  number  of 
the  flow,  Re  =  DU /v.  To  demonstrate  the  ability  of  the  code  to  resolve  the  laminar 
separation  at  low  Reynolds  numbers,  two-dimensional  solutions  were  computed  for 
Reynolds  numbers  between  15  and  40  in  increments  of  5. 

Solutions  were  computed  using  a  single  O-type  solution  domain  with  201  points 
in  the  circumferential  direction  and  201  points  in  the  radial  direction.  The  far  field 
boundary  was  located  at  95 D  from  the  cylinder  surface,  which  had  an  initial  spacing 
of  0.001.D.  Solutions  were  computed  using  the  HLLC  upwind  scheme  for  the  Euler 
terms  and  the  D3ADI  algorithm  with  Huang  sub-iteration  for  the  viscous  terms  with 
a  van  Albada  limiter.  Solutions  were  computed  using  a  freestream  Mach  number  of 
0.2  with  low  Mach  preconditioning  enabled  [61]. 

The  length  of  the  laminar  separation  bubble  produced  by  the  steady  simula¬ 
tions  was  determined  by  examining  the  u- velocity  profile  along  the  cylinder  centerline 
extending  in  the  freestream  direction  shown  in  Figure  8.  The  u- velocity  component 
is  negative  within  the  separation  bubble,  transitioning  to  positive  values  outside  the 
recirculation  region.  The  length  of  the  separation  bubble  expressed  in  units  of  L/D 
is  compared  to  experimental  values  and  previous  numerical  simulations  in  Figure  8. 
The  current  simulations  provide  good  agreement  to  the  previous  experimental  and 
numerical  result  presented  in  Bush  and  Baeder  [79]. 

3.1.2  Unsteady  Solution,  Re  >  35.  At  Reynolds  numbers  above  45,  flow 
instabilities  related  to  the  length  of  the  laminar  separation  bubble  lead  to  a  time- 
dependent  wake.  At  Reynolds  numbers  between  50  and  185,  the  stationary  circular 
cylinder  produces  a  two-dimensional  time-periodic  wake  structure  in  which  vortices 
of  alternating  strength  are  shed  from  the  cylinder  and  convect  downstream.  The 
shedding  frequency  of  this  von  Karman  vortex  street  can  be  expressed  using  a  non- 
dimensional  Strouhal  number,  St  =  fD/U,  which  exhibits  a  dependence  on  Reynolds 
number. 
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Figure  8:  L/D  as  a  function  of  Reynolds  number  (a)  and  the  u- velocity  profile 

behind  the  cylinder  as  a  function  of  x/D  (b).  Reference  data  digitized  from  Bush  and 
Baeder  [79]. 


Time-accurate  solutions  utilizing  the  same  single  block  solution  domain  de¬ 
scribed  for  the  steady  solutions  were  computed  at  Reynolds  numbers  60,  120,  and 
185  to  demonstrate  the  ability  of  the  low  Mach  number  preconditioning  to  capture 
a  time-dependent  flow  held.  Solutions  were  computed  using  a  non-dimensional  time 
step  of  0.005  with  20  dual  time  iterations  per  time  step. 

The  non-dimensional  shedding  frequencies  and  force  histories  predicted  by  OVER¬ 
FLOW  2.1-Elastic  exhibit  similar  trends  to  those  documented  in  previous  experimen¬ 
tal  and  numerical  investigations,  shown  in  Figure  9.  The  mean  drag  and  r.m.s.  lift  for 
the  solution  at  Re  185  are  compared  to  previous  numerical  solutions  of  Guilmineau  [80] 
presented  in  Table  1.  The  OVERFLOW  2.1-Elastic  solution  at  Reynolds  number  185 
over-predicts  mean  drag  and  r.m.s.  lift  by  8%  and  12%  when  compared  to  Lu  and 
Dalton  [81]. 


3.1.3  Prescribed  Transverse  Oscillation,  Re  =  185.  Dynamic  motion  of 
a  circular  cylinder  influences  the  vortex  wake  and  force  history  when  compared  to 
the  wake  produced  by  a  stationary  circular  cylinder.  An  oscillating  cylinder  will 
experience  a  phase  jump  in  the  vortex  wake  for  a  frequency  of  excitation,  fe,  near  the 
natural  shedding  frequency,  fQ  [80].  This  phenomenon  was  utilized  to  demonstrate 


28 


0.22 


0.2 

0.18 

0.16 

0.14 

0.12 

0.1 

40  60  80  100  120  140  160  180  200 


Williamson  (Experimental)  + 

Ponta  (Numerical)  X 

Baeder  (Numerical)  X 

Present  Solutions  □ 


□ 


#- 


4r 


£ 


Re 


(a) 


0  50  100  150  200  250  300  350  400 


CL  .  - 

■ 

0  50  100  150  200  250  300  350  400 


t*  (t  U/D) 

(b) 


Figure  9:  Strouhal  number  as  a  function  of  Reynolds  number  (a)  and  the  lift  and 
drag  history  (b).  Reference  data  digitized  from  Bush  and  Baeder  [79]. 


Table  1:  Comparison  of  the  mean  drag  and  r.m.s.  lift  produced  by  a  circular 

cylinder  at  Reynolds  number  185  to  previous  numerical  simulations. 

Source  Re  CDmean  CLrms  St  Method 


Present  Simulation 

60 

1.47 

0.088 

0.132 

Numerical 

Present  Simulation 

120 

1.41 

0.304 

0.171 

Numerical 

Present  Simulation 

185 

1.42 

0.477 

0.191 

Numerical 

Lu  and  Dalton  [81] 

185 

1.31 

0.422 

0.195 

Numerical 

Guilmineau  [80] 

185 

1.29 

0.443 

0.195 

Numerical 
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the  ability  of  OVERFLOW  2.1-Elastic  to  simulate  the  unsteady  flow  held  about  a 
moving  body. 

Solutions  were  computed  at  Reynolds  number  185  for  excitation  frequencies  of 
fe/fo  of  0.9  and  1.0  with  an  oscillation  amplitude,  Ae,  of  0.2 D  to  match  previous 
solutions  computed  by  Guilmincau  and  Queutey  [80].  The  behavior  of  the  lift  and 
drag  force  histories  match  those  computed  by  Guilmineau  and  Queutey  [80],  as  shown 
in  Figure  10  (b).  As  with  the  stationary  shedding  cylinder,  OVERFLOW  2.1-Elastic 
over  predicts  both  lift  and  drag.  However,  OVERFLOW  2.1-Elastic  does  capture 
the  trends  in  the  mean  drag  and  r.m.s  lift  exhibited  by  the  transversely  oscillating 
cylinder  at  different  frequencies,  shown  in  Figure  10  (a). 
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Figure  10:  Lift  and  drag  as  a  function  of  fe/f0  for  an  oscillation  amplitude  of  0.2 D, 
(a),  and  force  histories  (b).  Reference  data  digitized  from  Bush  and  Baeder  [79]. 


3.2  Two-Dimensional  Airfoil 

Additional  two-dimensional  dynamic  solutions  were  used  to  demonstrate  the 
ability  of  OVERFLOW  to  capture  the  time-dependent  how-held  generated  by  a  mov¬ 
ing,  or  dynamic  body  at  Reynolds  numbers  characteristic  of  the  three-dimensional 
happing  wing  problem.  Solutions  were  computed  to  match  experimental  data  for  a 
SD7003  airfoil  recorded  in  the  AFRL  Horizontal  Free-surface  Water  Tunnel.  Experi¬ 
mental  data  were  available  for  sinusoidal  heaving  motions  at  Reynolds  number  10,000. 
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The  data  set  includes  wake  visualization  using  dye  injection  at  the  leading  and  trail¬ 
ing  edges,  as  well  as  time-averaged  velocity  and  vorticity  measurements  obtained  with 
particle  image  velocimetry  (PIV)  [82].  Forces  were  not  measured  in  the  experimental 
investigation. 

The  two-dimensional  heaving  solutions  were  computed  for  a  constant  pitch  angle 
of  a{t)  =  4  degrees  and  a  non-dimensional  heaving  amplitude  of  h(t)  =  0.05  cos(u it)  at 
a  reduced  frequency  of  k  =  3.93.  Solutions  were  computed  at  a  chord-based  Reynolds 
number  of  10,000  to  match  the  experimental  conditions.  At  this  Reynolds  number  the 
flowfield  should  remain  laminar.  Solutions  were  computed  using  low  Mach  number 
preconditioning  at  a  freestream  Mach  number  of  0.015. 

The  experimental  data  set  does  not  include  force  measurements.  Thus,  the 
simulated  flowfield  was  qualitatively  compared  to  the  experimental  measurements 
of  velocity  magnitude  in  Figure  11  and  vorticity  magnitude  in  Figure  12.  The  time- 
dependent  wake  structures  predicted  by  the  OVERFLOW  2.1-Elastic  solution  achieve 
a  good  qualitative  match  to  the  wake  structures  documented  in  the  experimental  data 
set. 

The  heaving  airfoil  produced  small  vortical  structures  caused  by  flow  separation 
close  to  the  leading  edge  during  the  downward  phase  of  motion,  as  seen  in  Figure  13. 
As  the  wing  heaves,  the  vortex  convects  along  the  upper  surface,  slowing  about  the 
mid  chord,  then  accelerating  until  it  convects  into  the  freestream  from  the  trailing 
edge.  The  OVERFLOW  2.1-Elastic  solution  captures  this  unsteady  phenomenon, 
which  would  likely  be  missed  by  lower  order  models. 

3.3  3D  Flapping  Wing 

The  three-dimensional  configuration  is  based  on  the  geometry  and  dynamic  wing 
motion  of  the  Hawkmoth  at  hover.  Published  examinations  of  the  Hawkmoth  aero¬ 
dynamics  include  a  variety  of  numerical  and  experimental  investigations  that  have 
documented  the  characteristic  aerodynamic  phenomena.  However,  precise  experi- 
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Figure  11:  Contours  of  velocity  magnitude  normalized  by  freestream  velocity  for 

a  heaving  airfoil  predicted  using  OVERFLOW  2.1-Elastic,  left,  and  measured  using 
PIV  in  the  AFRL  Horizontal  Free-surface  Water  Tunnel,  right.  The  CFD  solution 
represents  a  single  snapshot  in  time  while  the  PIV  represents  a  mean  velocity  mag¬ 
nitude  over  multiple  snapshots  at  t/r  =  0.  The  non-dimensional  heaving  motion  is 
defined  as  h{t)  =  0.05cos(cat)  and  a(t)  =  4  for  a  reduced  frequency  of  k  —  3.93. 


32 


Figure  12:  Contours  of  vorticity  magnitude  normalized  by  freestream  velocity 

normal  to  the  cut  planes  for  a  heaving  airfoil  predicted  using  OVERFLOW  2.1-Elastic, 
left,  and  measured  using  PIV  in  the  AFRL  Horizontal  Free-surface  Water  Tunnel, 
right.  The  CFD  solution  represents  a  single  snapshot  in  time  while  the  PIV  represents 
a  mean  velocity  magnitude  over  multiple  snapshots  at  t/r  —  0.  The  non-dimensional 
heaving  motion  is  defined  as  h{t)  =  0.05cos(cuf)  and  a{t )  =  4  for  a  reduced  frequency 
of  k  =  3.93. 
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Figure  13:  Contours  of  vorticity  magnitude  normalized  by  freestream  velocity 

normal  to  the  cut  planes  for  a  heaving  airfoil  measured  using  PIV  in  the  AFRL 
Horizontal  Free-surface  Water  Tunnel  and  predicted  using  OVERFLOW  2.1-Elastic. 
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mental  measurements  of  the  aerodynamic  force  histories  produced  by  the  Hawkmoth 
during  flight  are  still  elusive.  Direct  measurement  of  the  aerodynamic  forces  are  com¬ 
plicated  by  the  small  force  amplitudes  and  high  frequencies,  as  well  as  the  need  to 
differentiate  the  inertial  elastic  forces  and  the  aerodynamic  forces  from  a  single  signal. 

The  lack  of  experimentally  determined,  time-resolved  aerodynamic  forces  limits 
the  validation  of  three-dimensional  flapping  wing  simulation  to  an  order  of  magnitude 
comparison  of  the  mean  value  of  the  simulated  force  history  and  the  force  required 
to  support  the  weight  of  a  typical  Hawkmoth  during  hover.  Beyond  the  single-point 
validation,  the  ability  of  the  Navier-Stokes  simulations  to  capture  the  documented  un¬ 
steady  aerodynamic  phenomena  provides  further  confidence  in  the  simulation  meth¬ 
ods. 

Numerical  simulations  may  also  be  compared  to  prior  simulation  efforts,  al¬ 
though  the  utility  of  a  comparison  to  prior  works  depends  on  the  specific  config¬ 
urations  being  compared.  Geometric  influences,  such  as  wing  planform  and  wing- 
body  configurations,  kinematic  variations,  and  assumptions  inherent  in  the  simulation 
method  can  all  influence  the  predicted  aerodynamic  forces.  Details  of  the  solution 
configuration  used  in  this  study  are  discussed  in  Chapter  IV. 

3.3.1  Test  Configuration.  The  characteristic  configuration  developed  for 
this  study  uses  a  symmetry  solution  for  a  single  wing  with  a  length  of  50  millimeters, 
offset  20  millimeters  from  the  centerline  symmetry  plane.  Two  wing  planforms  were 
considered,  a  Hawkmoth  wing  planform  based  on  the  work  of  Liu  [20],  and  a  rectan¬ 
gular  wing  with  NACA  0009  cross  section  and  a  30%  wing  chord,  shown  in  Figure  14. 
For  each  planform,  the  kinematic  rotations  are  prescribed  about  the  intersection  of 
the  leading  edge  and  the  wing  root.  The  large  centerline  offset  was  required  in  order 
to  resolve  the  rigid  body  rotations  while  preventing  the  trailing  edge  of  the  wing  from 
intersecting  the  idealized  body,  when  included,  or  the  symmetry  plane. 

The  rigid  wing  rotations  were  based  on  the  kinematic  parameters  documented 
by  Willmott  and  Ellington  [17,6]  for  a  tethered  Hawkmoth  at  hover.  These  kinematic 
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Table  2:  Current  and  published  mean  aerodynamic  forces  for  Hawkmoth  hover 

simulations.  Forces  are  presented  for  single  wing  symmetry  solutions.  Thrust  to 
weight  ratio,  T/W,  computed  based  on  a  mass  of  1.48  grams. 


Planform  Rectangular  Manduca  sexta  Agrius  convolvuli 


Study 

Current  Work 

Current  Work 

Aono  and  Liu  [22] 

Fx  (N) 

7.24xlO"04 

8.42xlO-04 

1.20xl0-03 

Fy  (N) 

-1.46xlO~03 

-1.65xlO-03 

-1.20x  10-03 

Fz  (N) 

7.26xl0-°3 

6.16xHT03 

8.48x  10~03 

T/W 

1.0 

0.8 

1.2 

parameters  formed  the  basis  for  previous  numerical  simulations  by  Liu  et  al.  [20, 19] 
and  Aono  et  al.  [22,24,25].  In  both  cases,  a  4-term  Fourier  approximation  of  the 
experimental  data  set  was  used  to  define  the  dynamic  wing  motions  for  the  numerical 
simulations. 

3.3.2  Aerodynamic  Forces.  Simulations  using  the  Liu  Hawkmoth  planform 
and  the  rectangular  planform  produced  mean  vertical  forces  able  to  sustain  the  mass 
of  a  Hawkmoth  at  hover,  as  outlined  in  Table  2.  These  mean  aerodynamic  forces 
are  comparable  to  the  published  simulation  of  Aono  and  Liu  [22],  However,  the  time 
history  of  the  vertical  forces  does  vary  between  the  simulations,  as  shown  in  Figure 
15.  The  time  resolved  aerodynamic  forces  computed  for  the  rectangular  wing  and 
the  Hawkmoth  planform  exhibit  a  somewhat  decreased  lift  amplitude  compared  to 
the  forces  computed  by  Aono  and  Liu  [22],  but  match  the  shape  of  the  force  history 
and  the  phase  of  the  peak  force  magnitudes.  In  particular,  the  current  simulations 
capture  the  “double  bump”  exhibited  during  both  half-strokes  of  the  4-term  rigid  wing 
simulations.  Variations  in  the  wing  planforms,  based  on  the  species  Agrius  convolvuli 
in  the  case  of  the  Aono  simulations,  the  location  of  the  wing  root,  or  simulation 
configuration  could  account  for  the  variation  in  the  force  magnitudes.  These  variations 
also  influence  the  time  history  of  the  aerodynamic  power,  shown  in  Figure  16. 

3.3.3  Wake  Structures.  The  ability  of  the  mean  vertical  force  component 
to  sustain  the  mass  of  the  Hawkmoth  at  hover  provides  a  single  point  validation 
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-  Agrius  convolvuli 

-  Manduca  sexta 

-  Rectangular 

Figure  14:  Wing  planforms  based  on  the  Manduca  sexta  used  by  Liu  and  Kawachi 
[20],  the  Agrius  convolvuli  used  by  Aono  and  Liu  [22],  and  the  Rectangular  wing 
planform  used  in  the  current  work. 
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Figure  15:  Time  resolved  rearward  and  vertical  force  histories  for  the  Hawkmoth 

at  hover.  Force  histories  represent  a  current  simulation  for  a  rectangular  planform,  a 
current  simulation  for  the  Manduca  sexta  planform,  and  digitized  forces  computed  by 
Aono  and  Liu  [22]  using  the  Agrius  convolvuli  planform. 
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Figure  16:  Time  resolved  aerodynamic  power  history  for  a  rectangular  planform,  a 
current  simulation  for  the  Manduca  sexta  planform,  and  digitized  aerodynamic  power 
computed  by  Aono  and  Liu  [22]  using  the  Agrius  convolvuli  planform. 

of  the  aerodynamic  forces  predicted  by  the  Navier-Stokes  simulation.  Comparing 
the  aerodynamic  phenomena  exhibited  by  the  time-accurate  Navier-Stokes  simulation 
to  previously  documented  aerodynamic  phenomenon  and  wake  structures  provides 
confidence  that  the  Navier-Stokes  simulation  is  capturing  the  relevant  aerodynamic 
phenomena. 

The  rectangular  and  Hawkmoth  planforms  each  produce  the  characteristic  wake 
structures  previously  documented  for  the  Hawkmoth  at  hover  using  experimental 
[15,16]  and  numerical  [19,24]  investigations.  The  leading  edge  vortex  and  tip  vortex 
documented  by  van  den  Berg  and  Ellington  [16]  are  shown  in  Figure  17  along  with 
vortical  wake  structures  produced  by  the  rectangular  wing  planform.  The  present 
simulations  capture  the  characteristic  horseshoe  wake  structure  formed  by  the  leading 
edge  vortex,  the  tip  vortex,  and  the  vortices  shed  at  the  start  of  the  half-stroke, 
highlighted  in  Figure  17.  The  smoke  visualizations  used  by  van  den  Berg  and  Ellington 
[15, 16]  highlight  the  leading  edge  vortex  and  the  evolution  of  the  tip  vortex,  but  do 
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not  highlight  the  other  wake  structures  such  as  the  starting  and  stopping  vortices  that 
may  be  seen  in  the  wake  structures  produced  by  the  present  simulations.  Examining 
the  wake  produced  by  the  rectangular  wing  planform  at  the  end  of  the  upstroke, 
shown  in  Figure  18,  shows  the  shed  circulation  at  the  trailing  edge  of  the  wing  as  well 
as  three-dimensional  influences  at  the  wing  root  produced  by  the  rigid  wing  rotations. 

3.3.4  Summary.  The  three-dimensional  simulations  based  on  the  rectangu¬ 
lar  wing  planform  and  the  Hawkmoth  planform  first  used  by  Liu  produce  the  aero¬ 
dynamic  forces  necessary  to  sustain  the  weight  of  a  Hawkmoth  at  hover.  In  addition, 
the  simulations  reproduce  the  characteristic  aerodynamic  phenomena  previously  doc¬ 
umented  for  the  Hawkmoth  at  hover. 
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Figure  17:  Leading  edge  vortex  and  tip  vortex  structures  documented  by  van  den  Berg  and  Ellington  [16]  using  a  scaled 
mechanical  flapping  mechanism,  shown  in  red,  and  the  wake  structures  produced  by  a  rigid  rectangular  planform,  below, 
using  contours  of  vortcity  magnitude.  The  vorticity  contours  show  all  the  near  body  wake  structures,  not  just  the  leading 
edge  vortex  and  tip  vortex  highlighted  using  smoke  visualization  by  van  den  Berg  and  Ellington. 


Figure  18:  Vortical  wake  structures  at  the  end  of  the  upstroke.  The  vortical  struc¬ 
tures  shed  at  the  beginning  of  the  stroke  can  be  seen  in  the  fore-field,  combining  with 
the  tip  vortex  and  the  leading  edge  vortex  to  form  a  horseshoe  vortex.  The  Leading 
edge  vortex  has  not  yet  shed  into  the  wake  as  the  wing  rotates  about  the  leading  edge. 
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IV.  Configuration  Study 

In  order  to  understand  how  geometric  variations  influence  the  predicted  aerodynamic 
performance  of  the  Hawkmoth  at  hover,  a  configuration  study  was  performed  that 
considered  the  influence  of  wing  planform,  wing-body  interactions,  and  wing  length 
on  the  time  resolved  aerodynamic  forces,  unsteady  aerodynamic  phenomena,  and  over¬ 
all  performance  metrics.  By  understanding  how  the  geometric  properties  influenced 
the  numerical  simulations,  a  baseline,  or  characteristic  configuration  was  developed 
which  reproduced  the  previously  documented  aerodynamic  attributes  for  the  Hawk- 
moth  at  hover  while  being  easily  reproducible  and  requiring  a  reasonable  amount  of 
computational  resources. 


4-1  Performance  Metrics 

In  order  to  understand  the  influence  of  configuration  on  performance,  simula¬ 
tions  were  examined  using  a  set  of  quantifiable  performance  metrics  based  on  the  com¬ 
puted  aerodynamic  forces  and  aerodynamic  power.  The  ability  of  the  configuration 
to  reproduce  the  characteristic  aerodynamic  phenomena  was  based  on  a  qualitative 
assessment  of  the  wake  structure  and  unsteady  aerodynamic  mechanisms. 


4-1.1  Aerodynamic  Forces  arid  Moments.  When  computing  time-accurate 
solutions,  OVERFLOW  2.1-Elastic  reports  aerodynamic  force  and  moment  coeffi¬ 
cients  as  a  function  of  time  step  in  the  fomoco.out  output  file  for  each  integration 
surface  defined  in  the  MIXSUR  input  file,  mixsur.inp.  The  force  and  moment  coef¬ 
ficients  reported  in  the  FOMOCO  output  are  non-dimensionalized  by  the  reference 
velocity  defined  in  the  OVERFLOW  input  deck  and  the  surface  area  defined  for  each 
integration  surface  in  the  MIXSUR  input. 
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For  dynamic  solutions,  the  position,  orientation,  and  non-dimensional  velocity 
of  each  moving  component  are  recorded  in  the  ANIMATE  output  hie,  animate,  out. 
The  ANIMATE  output  also  includes  the  non-dimensional  forces  and  moments  for 
each  moving  component. 

F 

Fan  i  mate  =  ~  772  72  (8) 

PrefU  reflref 


M ANIMATE  ~  tt2  t  3  (9) 

PrefU  reflref 

From  this  data,  mean  and  peak  aerodynamic  forces  may  be  computed.  In  addi¬ 
tion,  the  time  resolved  forces  may  be  combined  with  the  prescribed  wing  kinematics  to 
examine  the  phase  relationships  between  the  aerodynamic  forces  and  the  prescribed 
wing  motions. 


4-1.2  Aerodynamic  Power  and  Efficiency.  For  dynamic  solutions,  the  AN¬ 
IMATE  output  may  be  used  to  compute  the  aerodynamic  power  for  each  moving 
component.  The  aerodynamic  power  required  to  actuate  a  rigid  body  may  be  com¬ 
puted  using  the  inner  product  of  the  aerodynamic  forces  and  moments  with  the  body’s 
translational  and  rotational  velocities  using  Equation  (10). 
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The  aerodynamic  power  may  be  used  to  compute  an  aerodynamic  efficiency  based  on 
the  ratio  of  vertical  aerodynamic  force  component  to  aerodynamic  power. 
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This  efficiency  is  analogous  to  the  definition  of  Power  Loading,  PL,  used  in  the  rotor 
craft  community  [83]. 
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4-1.3  Inertial  Power.  Included  in  the  ANIMATE  output  are  the  rotational 
velocities  of  each  moving  component  in  the  body  reference  frame.  This  allows  for 
a  simplified  computation  of  inertial  power  requirements  for  rigid  body  motion.  The 
inertial  power  for  the  rigid  body  wing  motion  may  be  computed  using  the  inertial 
moment  at  the  wing  root  and  the  rotational  velocities  due  to  the  wing  motion. 


Pinertial  A/ inertial  '  &wing  T  Finertial  '  V wing 


(12) 


The  moment  at  the  wing  root  due  to  inertial  forces  for  a  rigid  wing  may  be  computed 
using  the  inertia  tensor  and  rotational  accelerations  due  to  body  motion. 
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The  inertia  tensor,  [I],  represents  the  mass  distribution  of  the  specific  wing  configu¬ 
ration  being  studied.  The  center  of  mass  of  various  insect  wings  have  been  previously 
documented  [12,30],  which  may  be  used  to  approximate  an  inertia  tensor.  In  the  case 
of  a  flexible  wing  structure,  the  inertia  tensor  will  be  time  variant  with  the  structural 
response.  As  with  the  aerodynamic  power,  the  inertial  power  may  be  used  to  compute 
an  inertial  efficiency  metric. 

_\Finertial\ 
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4-1-4  Aerodynamic  Force  and  Moment  as  a  Function  of  Span.  Included  as 
part  of  the  fluid-structure  coupling  framework  in  OVERFLOW  2.1-Elastic,  aerody¬ 
namic  forces  and  moments  are  reported  as  a  function  of  span  along  the  clastic  axis 
defined  for  each  flexible  body.  These  aerodynamic  forces  may  be  reported  as  lumped 
aerodynamic  forces  for  each  wing  section  or  as  a  force  per  unit  span  based  on  the 
&EROTOR  namelist  inputs.  The  ROTOR  forces  are  non-dimensionalized  by  ref- 
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erence  parameters  specified  in  the  &EROTOR  namelist  in  the  OVERFLOW  input 
deck. 


ROTOR  g  _  5pROTOR  (jl  rot  OrM rotor  )2LreF 
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Visualizing  the  aerodynamic  forces  as  a  function  of  span  provides  insight  into 
the  force  distribution  along  the  wing  surface,  and  provides  insight  into  which  regions 
of  the  wing  contribute  to  the  vertical  and  side  force  production  throughout  the  wing 
stroke. 


4-1.5  N on- dimensionalization  Based  on  Tip  Velocity.  Non-dimensional 
aerodynamic  forces  and  moments  as  defined  above  are  based  on  an  arbitrary  ref¬ 
erence  velocity  which  does  not  have  significant  meaning  as  the  freestream  velocity 
approaches  zero.  However,  an  alternate  non-dimensionalization  for  a  flier  in  a  hover 
condition  may  be  formulated  based  on  the  maximum  tip  velocity  of  the  oscillating 
wing. 
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This  non-dimensionalization  has  a  corollary  for  rotorcraft,  where  it  is  referred  to  as 
an  area  weighted  thrust  coefficient  [83].  Unless  otherwise  specified,  the  reader  may 
assume  that  the  tip  velocity  used  in  the  non-dimensionalization  is  based  on  a  sinusoidal 
approximation  of  the  wing  kinematics,  given  in  Equation  (20),  based  on  the  flapping 
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frequency,  /,  peak-to-peak  sweep  amplitude,  $,  and  wing  length,  R. 


Vtip  =  (2vr/)  (0.5$)  R 


(20) 


4-2  Wing  Planform 

The  two  planforms  considered  for  this  study  were  the  Hawkmoth  wing  planform 
based  on  wing  geometry  published  by  Liu  [20],  and  an  idealized  rectangular  wing  with 
a  30%  chord  based  on  wing  length.  The  surface  domains  used  to  model  the  left  wing 
for  each  planform  are  shown  in  Figure  19.  The  simulations  presented  were  computed 
for  the  same  torsionally  flexible  configuration  with  the  wing  root  and  pivot  location 
offset  from  the  centerline  by  40%  of  the  wing  length  described  in  Section  4.3.  The 
two  planform  configurations  with  wing  motion  are  shown  in  Figure  20. 

The  Hawkmoth  planform  defined  by  Liu  represents  an  approximation  of  the 
two  wing  pairs  that  overlap  during  the  wing  stroke  to  form  a  single  aerodynamic 
surface.  During  the  wing  stroke,  the  two  wing  pairs  exhibit  small  relative  movements 
that  influence  the  shape  of  the  wing  planform,  the  area  of  the  combined  aerodynamic 
surface,  and  the  relative  incidence  of  the  two  wings.  Capturing  the  dynamic  behavior 
into  a  single  wing  represents  a  best  approximation  of  this  dynamic  wing  structure. 
The  Liu  planform  represents  an  approximation  of  a  particular  wing  pair,  which  varies 
with  time  and  exhibits  small  variations  between  individual  Hawkmoth  samples.  The 
rectangular  wing  represents  a  simplified  planform  that  is  easily  reproducible  for  future 
numerical  and  experimental  investigations. 

The  time-resolved  aerodynamic  forces  for  the  rectangular  planform  and  the 
Hawkmoth  wing  planform  are  presented  in  Figure  21.  The  rectangular  wing  and 
the  Hawkmoth  wing  planforms  both  produce  force  histories  of  a  similar  character, 
but  with  deviations  of  the  peak  vertical  forces  during  each  half  stroke  of  up  to  30%. 

The  vortical  wake  structures  for  the  rectangular  wing  planform  and  the  Hawk¬ 
moth  wing  planform  are  illustrated  in  Figure  22.  Both  wing  planforms  produce  similar 
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Figure  19:  Comparison  of  Hawkmoth  wing  planform,  upper,  with  the  rectangular 
wing  planform,  lower.  Wings  correspond  to  the  left  wing  when  viewed  from  above. 
Each  wing  consists  of  three  overset  solution  domains.  For  the  Hawkmoth  planform, 
the  domains  are  a  collar  (blue),  an  upper  skin  (red),  and  a  lower  skin  (not  shown). 
For  the  rectangular  wing,  the  solution  domains  are  a  main  wing  skin,  shown  in  red, 
an  outer  wing  cap,  and  an  inner  wing  cap,  each  shown  in  blue. 
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Figure  20:  Illustration  of  the  Hawkmoth  wing  planform  configuration,  (a)  and 

(b),  and  the  rectangular  wing  planform  configuration,  (c)  and  (d).  The  direction  of 
the  wing  motion  is  indicated  using  blue  arrows.  Wing  motions  represent  the  4-term 
kinematic  approximation  with  prescribed  feathering  along  the  wing  length. 
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overall  wake  structures,  reproducing  the  character  of  the  unsteady  Hawkmoth  wake. 
However  there  are  local  small  variations  in  the  magnitude  and  the  strength  of  the 
wake  structures  produced  by  the  two  planforms. 

4-3  Wing-Body  Interactions 

The  Hawkmoth  wing  is  a  flexible  structure  that  is  joined  to  the  body  with 
a  compliant  membrane  at  the  wing  root.  It  is  not  possible  to  model  this  flexible 
membrane  and  preserve  the  observed  wing-body  interactions  when  using  a  rigid  wing 
simulation.  The  ability  to  examine  the  performance  of  a  rigid  wing  necessitated  that 
the  wing  root  be  offset  from  the  body  in  order  to  resolve  the  motion  of  the  wing  root 
while  preventing  contact  between  the  wing  surface  and  the  body  surface. 

To  understand  the  influence  of  wing-body  interactions  on  the  unsteady  aerody¬ 
namic  mechanisms,  simulations  were  computed  with  and  without  an  idealized  Hawk¬ 
moth  body.  The  Hawkmoth  body  was  represented  using  an  elliptical  body  with  a 
circular  cross  section.  The  body  length  was  80%  of  the  wing  length  and  the  body 
diameter  was  10%  of  the  wing  length.  The  simulations  presented  were  computed  for 
torsionally  flexible  wings  with  the  wing  root  and  pivot  location  offset  from  the  cen¬ 
terline  by  40%  of  the  wing  length.  The  40%  centerline  offset  was  the  minimum  offset 
that  allowed  full  kinematic  motion  of  the  rigid  Hawkmoth  wing  while  maintaining  the 
grid  spacing  requirements  in  the  regions  where  the  wing  and  body  solution  domains 
overlapped  in  physical  space.  For  wing  offsets  less  than  40%,  the  wing  root  trailing 
edge  of  the  rigid  Hawkmoth  wing  would  pierce  the  body  surface  at  pronation  and 
supination,  which  correspond  to  the  maximum  wing  stroke  angles.  Smaller  wing  off¬ 
sets  may  be  achieved  for  alternative  wing  planforms  or  for  flexible  wings,  depending 
on  the  specific  planform  or  set  of  wing  deformations.  The  solution  domains  used  to 
compute  the  flexible  rectangular  wing-body  configuration  arc  illustrated  in  Figure  23. 
The  solution  configuration  for  the  wing-only  simulations  utilized  the  same  near-body 
and  off-body  solution  domains  with  the  exception  of  the  body  domain. 
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Figure  21:  Comparison  of  the  rectangular  wing  solution  with  elliptical  body  to 

the  Hawkmoth  planform  simulations  with  elliptical  body.  Simulations  are  based  on  a 
4-term  kinematic  approximation  with  prescribed  feathering  along  the  wing  length. 
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Figure  22:  Comparison  of  the  wake  structures  of  the  rectangular  wing  planform 

to  the  Hawkmoth  wing  planform.  Simulations  are  based  on  a  4-term  kinematic  ap¬ 
proximation  with  prescribed  feathering  along  the  wing  length.  Wake  structures  are 
iso-surface  of  vorticity  magnitude. 
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Figure  23:  Near-body  solution  domains  (a)  and  off-body  solution  domains  (b)  for 

the  flexible  wing-body  configuration  at  t/r  =  0.25.  The  solution  domains  representing 
the  wing  are  colored  in  blue  and  green  while  the  solution  domain  representing  the 
elliptical  body  is  colored  in  red.  In  (a),  the  first  off-body  domain  is  represented  in 
grey,  while  the  set  of  off  body  domains  are  represented  in  black  in  (b). 
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The  computed  aerodynamic  forces  for  the  two  rectangular  wing  simulations, 
with  and  without  the  elliptical  body,  are  shown  in  Figure  24.  The  force  histories 
are  almost  identical,  indicating  that  the  body  does  not  significantly  impact  the  aero¬ 
dynamic  mechanisms  generated  by  the  wing  motions  for  simulations  with  the  large 
wing-body  gap.  The  wake  structures  produced  by  the  two  configurations  show  mini¬ 
mal  variation  between  the  wing-body  and  the  wing  only  simulations,  shown  in  Figure 
25. 


4-4  Influence  of  Wing  Length 

To  understand  the  influence  of  wing  length  on  the  unsteady  aerodynamic  mech¬ 
anisms,  flow  solutions  were  computed  for  three  wing  lengths  utilizing  identical  wing 
kinematics.  Solutions  were  computed  for  dimensional  wing  lengths  of  30,  50,  and  80 
mm,  with  a  constant  aspect  ratio.  It  was  observed  that  the  magnitude  of  the  di¬ 
mensional  forces  and  moments  produced  by  the  three  wing  lengths,  shown  in  Figure 
26,  scaled  with  the  square  of  the  tip  velocity,  Vfp.  When  the  forces  and  moments 
are  non-dimensionalized  based  on  the  tip  velocity,  using  Equations  (17)  and  (18), 
the  non-dimensional  force  histories  reduce  to  similar  time  histories  of  same  order  of 
magnitude,  as  shown  in  Figure  27.  The  non-dimensional  force  histories  do  show  small 
variations  in  the  peak  lift  amplitude.  It  is  possible  that  the  force  variations  are  due  to 
variations  in  the  Reynolds  numbers  used  for  each  simulation.  The  Reynolds  number 
used  for  the  three  simulations  was  held  constant,  and  not  updated  for  the  change 
in  tip  velocity  due  to  the  change  in  wing  length.  That  the  time  resolved  aerody¬ 
namic  forces  produced  by  wings  of  different  lengths  may  be  non-dimensionalized  by 
tip  velocity  is  significant  in  that  it  removes  wing  length  from  the  design  space  to 
be  considered  when  examining  the  influence  of  kinematic  variations  and  structural 
response  on  aerodynamic  performance. 
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Figure  24:  Comparison  of  the  wing-body  solution  to  the  wing-only  solution  for  the 

rectangular  planform.  Simulations  are  based  on  a  4-term  kinematic  approximation 
with  prescribed  feathering  along  the  wing  length.  The  time-resolved  aerodynamic 
force  coefficients  are  not  significantly  influenced  by  the  presence  of  the  elliptical  body. 
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Figure  25:  Comparison  of  the  wake  structures  of  the  rectangular  wing  solution 

with  and  without  the  elliptical  body.  Simulations  are  based  on  a  4-term  kinematic 
approximation  with  prescribed  feathering  along  the  wing  length.  Wake  structures  are 
iso-surface  of  vorticity  magnitude. 
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Figure  26:  Comparison  of  the  dimensional  forces  computed  for  the  30  mm,  50  mm, 
and  80  mm  rectangular  wings  with  prescribed  torsion. 
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Figure  27:  Comparison  of  the  force  coefficients  based  on  tip  velocity  computed  for 
the  30  mm,  50  mm,  and  80  mm  rectangular  wings  with  prescribed  torsion. 


4 . 5  Summary 

The  aerodynamic  forces,  moments,  and  power  requirements  non-dimensionalized 
based  on  the  approximate  peak  tip  velocity  provide  a  set  of  performance  metrics  that 
may  be  used  to  compare  solutions  for  varying  configurations. 

Of  the  geometric  parameters  considered,  wing-body  interaction,  wing  planform, 
and  wing  length,  wing  planform  had  the  most  significant  influence  on  the  unsteady 
aerodynamic  mechanisms.  The  influence  of  wing-body  interaction  was  minimal  for 
the  size  of  the  gap  necessary  to  resolve  the  rigid  wing  motion  without  contacting  the 
body  surface.  It  was  shown  that  the  influence  of  wing  length  on  the  aerodynamic 
performance  could  be  accounted  for  by  non-dimensionalizing  the  aerodynamic  forces 
by  the  square  of  the  peak  of  the  wing  tip  velocity. 

The  influence  of  the  wing  planform  could  be  seen  in  the  time-resolved  aerody¬ 
namic  forces  and  vortical  wake  structures  produced  during  the  wing  stroke.  Unfor¬ 
tunately,  the  Hawkmoth  planform  used  in  this  study  is  a  rough  approximation  of  an 
actual  Hawkmoth  wing,  and  does  not  provide  a  good  representation  of  the  wing-body 
junction.  Realizing  that  the  Hawkmoth  planform  was  a  rough  approximation  of  an 
actual  Hawkmoth  wing  and  that  the  rectangular  wing  planform  captures  the  relevant 
unsteady  aerodynamic  mechanisms  exhibited  by  the  Hawkmoth  in  flight,  the  baseline 
configuration  used  for  the  remainder  of  the  study  is  based  on  a  symmetry  solution 
with  a  single  rectangular  wing  planform  and  no  body,  as  described  in  Section  3.3. 
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V.  Analysis  of  Kinematic  Variations 

Due  to  the  correlation  between  the  unsteady  aerodynamic  phenomena  and  the  dy¬ 
namic  wing  motion,  variations  in  the  dynamic  wing  motions  have  a  direct  influence  on 
the  unsteady  aerodynamic  mechanisms  and  aerodynamic  performance  [2, 3, 4, 5, 6,  7]. 
Thus,  it  is  important  to  have  a  basic  understanding  of  the  influence  of  the  kinematic 
profile  on  the  unsteady  aerodynamic  phenomena  before  examining  the  influence  of 
structural  response  on  aerodynamic  performance. 

The  motion  of  a  rigid  flapping  wing  may  be  defined  using  a  set  of  time-dependent 
translations  and  rotations  about  the  wing  root  using  any  of  the  standard  engineering 
notations.  There  are,  however,  alternative  notations  that  describe  the  wing  motion 
based  on  the  physiology  of  a  given  flier,  often  referred  to  in  literature  as  kinematic 
parameters.  The  choice  of  kinematic  parameters  used  to  describe  the  wing  motion 
of  a  given  organism  may  be  based  on  skeletal,  musculature,  or  other  physiological 
constraints  of  the  organism,  or  they  may  be  defined  for  convenience  in  the  collection 
of  experimental  data. 

5.1  The  Kinematic  Parameters 

The  kinematic  parameters  used  for  this  study  to  define  the  rigid  wing  motion 
are  based  on  the  experimental  work  of  Willmott  and  Ellington  [84,17,6].  Willmott 
and  Ellington  documented  the  time-dependent  motion  of  the  Hawkmoth  wing  us¬ 
ing  four  kinematic  rotations.  The  four  kinematic  parameters  that  define  the  wing 
orientation  in  a  global  reference  frame  are  sweep,  </>(£),  elevation,  9(t),  feathering, 
a(t ),  and  stroke-plane  inclination,  /3,  as  shown  in  Figure  28.  The  stroke-plane  passes 
through  the  wing  root  and  is  inclined  with  respect  to  the  freestream  orientation.  The 
stroke-plane  represents  a  reference  plane  that  captures  the  dominant  sweep  rotation, 
as  shown  in  Figure  29.  The  sweep  angle,  represents  the  angle  of  the  wing  axis, 
or  feathering  axis,  projected  onto  the  stroke-plane.  Wing  sweep  is  often  referred  to 
as  the  in-plane  rotation.  The  elevation  angle,  9{t),  represents  the  deviation  from  the 
stroke-plane,  or  the  angle  measured  between  the  feathering  axis  and  its  projection 
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onto  the  stroke-plane.  The  elevation  angle  is  also  known  as  the  out-of-plane  rotation. 
Together,  the  sweep  and  elevation  angle  define  the  three-dimensional  motion  of  the 
wing  tip.  The  feathering  angle,  a(t),  represents  the  rotation  about  the  feathering  axis. 
The  feathering  axis  corresponds  to  the  leading  edge  of  the  wing  for  the  rectangular 
wing  simulations,  and  a  line  from  the  wing  root  to  the  wing  tip  for  simulations  using 
the  Hawkmoth  planform.  A  fifth  angle,  defines  the  orientation  of  the  body  with 
respect  to  the  freestream  direction  in  the  inertial  reference  frame.  Body  orientation 
as  defined  does  not  directly  influence  the  dynamic  wing  motion,  but  it  does  influence 
the  geometric  configuration  and  wing-body  relationship.  For  this  analysis,  sweep, 
elevation,  and  feathering  angles  vary  as  a  function  of  time,  while  stroke-plane  incli¬ 
nation  and  body  angle  are  constants.  Together,  these  five  parameters  describe  the 
orientation  of  the  body  and  the  rigid  body  motion  of  the  wing  in  the  inertial  reference 
as  a  function  of  time. 


Figure  28:  Hawkmoth  kinematic  parameters  defined  by  Willmott  and  Ellington  [17]. 
The  sweep  angle,  0(f),  represents  the  orientation  of  the  wing  projected  onto  the  stroke- 
plane.  The  elevation  angle,  0(f),  represents  the  angle  of  the  wing  measured  between 
the  feathering  axis  and  the  stroke-plane,  and  the  feathering  angle,  a(f),  represents 
the  rotation  about  the  feathering  axis. 
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Figure  29:  The  leading  edge  tip  path  of  the  left  wing,  shown  in  black,  with  respect 

to  the  stroke-plane,  shown  in  green.  The  wing  orientation  during  the  stroke  cycle  is 
shown  in  red  with  blue  circles  indicating  the  leading  edge. 

5.1.1  Phases  of  Motion.  When  examining  the  influence  of  the  kinematic 
profile,  it  is  convenient  to  decompose  the  wing  motion  into  phases  based  on  the  type 
of  motion.  The  wing  motion  during  the  flapping  cycle  is  composed  of  two  dominant 
types  of  motion:  a  wing  sweep  in  which  the  wing  moves  “up”  or  “down”,  and  a 
feathering  rotation  about  the  longitudinal,  or  span-wise,  wing  axis.  The  use  of  “up” 
and  “down”  to  define  the  direction  of  wing  motion  is  often  misleading,  as  “up”  in 
this  case  refers  to  the  upstroke  based  on  the  Hawkmoth  morphology.  The  direction  of 
the  upstroke  in  the  inertial  reference  frame  depends  on  the  inclination  of  the  stroke- 
plane,  which  is  only  15  degrees  from  freestream  direction  for  the  Hawkmoth  at  hover, 
as  shown  in  Figure  29.  The  wing  sweep  is  characterized  by  large  tip  velocities  as  the 
wing  sweeps  through  111  degrees.  The  wing  sweep,  or  “translational  motion”,  and 
the  feathering,  or  rotational  motion,  are  not  separated  in  time;  rather  they  overlap  at 
the  end  of  one  half-stroke  and  the  beginning  of  the  following  half-stroke  as  shown  in 
Figure  30.  For  the  Hawkmoth,  the  wing  feathers  through  120  degrees  at  the  end  of 
each  sweep  cycle  such  that  the  same  edge  of  the  wing  leads  during  both  the  up-  and 
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downstrokes.  With  a  flapping  frequency  of  26  Hz,  the  wing  motion  is  dominated  by 
large  rotational  rates  and  accelerations  during  each  flapping  cycle. 

5.1.2  Kinematic  Parameters  to  Rigid  Wing  Motion.  In  order  to  utilize  the 
kinematic  parameters  to  prescribe  the  wing  motion  for  the  numerical  simulations,  a 
direction  cosine  matrix  was  computed  by  applying  the  kinematic  parameters  in  the 
following  order: 

1.  Apply  the  feathering  rotation,  cc,  about  wing-local  x-axis 

2.  Apply  the  elevation  rotation,  9 ,  about  body-fixed  y- axis 

3.  Apply  the  sweep  rotation,  0,  about  body-fixed  z-axis 

4.  Apply  the  stroke-plane  inclination  rotation  ,  f3,  about  body-fixed  x-axis 

The  direction  cosine  matrix  may  then  be  used  to  transform  the  wing  from  an  initial 
location  to  the  time-dependent  orientation,  or  the  direction  cosine  matrix  may  be 
used  to  compute  the  quaternions  utilized  by  OVERFLOW  2.1-Elastic  to  define  wing 
orientation  as  a  function  of  time.  This  approach  was  used  to  compute  the  wing  orien¬ 
tation  throughout  the  wing  stroke  cycle  illustrated  in  Figure  30  (a)  for  the  kinematic 
parameters  in  Figure  30  (b). 

5.2  Influence  of  the  Kinematic  Approximations  on  Aerodynamic  Force 
Production 

In  order  to  remove  discontinuities  in  the  experimentally  determined  kinematic 
parameters  and  interpolate  those  parameters  at  a  time  resolution  suitable  for  numeri¬ 
cal  simulation,  a  truncated  Fourier  series,  as  defined  in  Equation  (21),  was  utilized  to 
approximate  each  of  the  three  time-dependent  kinematic  parameters  given  in  Equa¬ 
tions  (22),  (23),  and  (24). 


(21) 


n=l 
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Figure  30:  Rigid  wing  motion  in  the  inertial  reference  frame  (a)  for  the  Hawkmoth 
hover  kinematics  used  by  Liu  [19,20]  and  Aono  [24]  (b).  The  outline  of  the  Hawkmoth 
body  is  included  for  orientation. 
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Both  Liu  and  Aono  utilized  a  four-term  Fourier  approximation  to  smooth  the  exper¬ 
imental  kinematic  data  measured  by  Willmott  and  Ellington.  However,  the  influence 
of  the  truncated  Fourier  approximation  on  the  aerodynamic  performance  of  the  re¬ 
sulting  numerical  simulations  was  not  addressed  in  the  available  literature. 
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Prior  investigations  examining  the  influence  of  kinematic  variations  on  rigid 
wing  performance  by  Bos  [85,  86,  87,  88]  indicate  that  wing  kinematics  that  mimic 
the  observed  kinematics  of  insects  in  flight  produce  greater  lift-to-drag  ratios  than 
simplified  or  harmonic  kinematic  approximations.  This  has  a  direct  influence  on 
aerodynamic  power  requirements.  Other  investigations  examining  the  influence  of 
wing  kinematics  [89,90,91,7,2,3,4,5,6]  support  similar  conclusions  that  simplified 
wing  kinematics  are  often  unable  to  recreate  the  aerodynamic  performance  of  natural 
flapping  wing  fliers.  Thus,  it  is  important  to  understand  the  influence  of  the  kine¬ 
matic  approximation  on  aerodynamic  performance  prior  to  examining  the  influence 
of  structural  flexibility  on  aerodynamic  performance. 

A  full  sensitivity  study  examining  the  influence  of  kinematic  variation  is  outside 
the  scope  of  the  current  research.  However,  by  varying  the  number  of  terms  used  in 
the  Fourier  approximation  of  the  kinematic  parameters,  the  sensitivity  of  the  aerody¬ 
namic  mechanisms  and  aerodynamic  performance  to  the  Fourier  approximation  may 
be  examined,  and  by  extension,  lead  to  a  greater  understanding  of  the  role  of  the 
kinematic  parameters  and  kinematic  rates  in  determining  aerodynamic  performance. 
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5.2.1  The  Fourier  Approximation  of  the  Kinematic  Parameters.  To  assess 
the  influence  of  the  kinematic  approximation  on  aerodynamic  performance,  simula¬ 
tions  were  computed  for  hover  wing  kinematics  published  by  Aono  [24]  approximated 
using  2-,  3-,  and  4-term  Fourier  series.  The  4-term  Fourier  series  represents  a  complete 
reconstruction  of  the  kinematic  parameters  used  by  Liu  and  Aono,  as  their  kinematic 
parameters  were  based  on  a  4-term  approximation  of  the  experimentally  measured 
kinematic  parameters  of  Willmott  and  Ellington  [17,6].  The  approximated  kinematic 
parameters  and  kinematic  rates  for  the  2-,  3-,  and  4-term  Fourier  approximations  are 
shown  in  Figure  31  along  with  the  digitized  points  representing  the  4-term  kinematic 
profile  used  by  Aono. 

Comparing  the  2-term  Fourier  approximation,  or  harmonic  approximation,  of 
the  kinematic  parameters  in  Figure  31  to  the  digitized  4-term  approximation,  the  2- 
term  approximation  captures  the  amplitude  of  the  flapping  and  feathering  rotations 
but  does  not  capture  the  out-of-plane  rotation.  The  2-term  kinematic  approximation 
also  preserves  the  phase  relationship  between  the  sweep  and  feathering  parameters, 
but  does  not  capture  the  phase  relationship  between  the  peak  sweep  and  feathering 
rates.  It  is  not  until  all  four  terms  are  included  in  the  Fourier  approximation  that 
the  phase  relationship  between  the  peak  sweep  and  feathering  rates  used  by  Aono  are 
captured. 

Because  the  2-term  approximation  represents  a  harmonic  approximation  of  the 
kinematic  parameters,  it  is  not  able  to  capture  the  shortened  duration  of  the  upstroke 
compared  to  the  downstroke,  nor  does  it  capture  the  squared-off  shape  of  the  feather¬ 
ing  curve.  Compared  to  the  2-term  approximation,  the  3-term  approximation  captures 
a  bias  in  the  duration  of  the  upstroke  to  downstroke,  and  provides  a  slightly  better 
capture  of  the  sweep  rate.  The  3-term  approximation  also  captures  the  out-of-plane 
rotation  as  well  as  the  elevation  rate.  However,  the  3-term  approximation  still  does 
not  capture  the  squared  off  shape  of  the  feathering  parameter  or  the  complex  shape 
of  the  feathering  rate.  The  4-term  approximation  provides  a  complete  reconstruc¬ 
tion  of  the  kinematic  parameters  as  used  by  Liu  and  Aono.  Compared  to  the  2-  and 
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3-term  Fourier  approximations,  the  4-term  approximation  exhibits  a  complex  phase 
relationship  between  the  sweep,  elevation,  and  feathering  rates  which  influences  the 
tip  velocity  of  the  wing  and  the  duration  of  the  feathering  rotation  during  pronation 
and  supination. 

5.2.2  Influence  of  Kinematic  Variations  on  Aerodynamic  Performance.  The 
rigid  wing  simulations  based  on  the  2-,  3-,  and  4-term  Fourier  approximations  of  the 
wing  kinematics  each  produce  distinct  aerodynamic  force  histories,  shown  in  Figure 
32.  The  2-  and  3-term  rigid  wing  solutions  both  produce  force  histories  with  a  single 
peak  in  the  force  magnitude  late  in  each  half-stroke,  while  the  4-term  rigid  wing 
solution  produces  a  force  history  with  two  peaks  per  half-stroke.  Though  the  rigid 
wing  solutions  for  the  2-  and  3-term  kinematic  profiles  both  produce  a  single  peak  in 
the  aerodynamic  force  magnitude,  the  amplitude  and  history  of  the  directional  forces 
are  unique  to  each  kinematic  profile.  For  example,  the  2-term  rigid  wing  solution  is 
the  only  solution  presented  that  generates  a  negative,  or  downwards  vertical  force  at 
any  point  in  the  flapping  cycle.  Yet,  the  2-term  kinematic  profile  also  generates  peak 
vertical  forces  that  are  almost  twice  the  amplitude  of  the  vertical  forces  generated  by 
the  3-term  kinematic  profile. 

To  understand  why  the  three  kinematic  approximations  produce  the  distinct 
force  histories,  the  relationship  between  the  aerodynamic  forces  and  the  kinematic  pa¬ 
rameters  was  examined  by  plotting  the  aerodynamic  forces  against  the  time-dependent 
kinematic  parameters,  kinematic  rates,  and  kinematic  accelerations,  as  shown  in  Fig¬ 
ures  33,  34,  and  35.  These  phase  plots  highlight  correlations  between  the  kinematic 
parameters  and  the  aerodynamic  forces. 

Examining  the  phase  plots,  the  relationship  between  the  peak  aerodynamic  force 
magnitude  and  the  sweep  and  elevation  rates  stand  out  in  Figure  34.  Based  on  an 
understanding  of  unsteady  aerodynamic  mechanisms  and  delayed  stall,  it  is  a  logical 
assumption  that  the  wing  should  produce  peak  aerodynamic  force  during  the  trans¬ 
lation  phase  when  the  wing  reaches  peak  translational  velocity  at  a  moderate  angle 
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Figure  31:  Smoothed  kinematic  parameters  and  kinematic  rates  based  on  2-term,  3- 
term,  and  4-term  Fourier  approximations  of  the  published  Aono  kinematic  parameters 
[24].  The  digitized  Aono  parameters  were  based  on  a  4-term  Fouier  approximation  of 
the  Willmott  and  Ellington  parameters  [17]. 
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of  attack.  Examining  the  phase  relationship  between  aerodynamic  force  magnitude 
and  the  sweep  rate  for  the  2-,  3-,  and  4-term  kinematic  profiles  shown  in  Figure  34, 
the  correlation  between  peak  force  magnitude  and  sweep  rate  was  strongest  for  the 
4-term  kinematic  profile,  while  the  peak  aerodynamic  force  magnitude  lagged  behind 
the  peak  sweep  rate  for  the  2-  and  3-term  kinematic  profiles.  In  Figure  34,  the  kine¬ 
matic  rates  are  plotted  along  the  horizontal  axis  and  the  force  magnitude  is  plotted  on 
the  vertical  axis.  When  the  peak  aerodynamic  force  magnitude  occurs  in  phase  with  a 
peak  kinematic  rate,  the  plot  of  the  force  history  forms  a  sharp  peak,  as  exhibited  by 
the  4-term  profile.  When  the  peak  kinematic  rate  precedes  the  peak  force  magnitude, 
the  force  plot  forms  an  open  loop  as  exhibited  by  the  3-term  kinematic  profile. 

The  3-  and  4-term  kinematic  profiles  both  exhibited  strong  correlations  between 
peak  aerodynamic  forces  and  peak  elevation  rate,  while  the  2-term  profile  exhibited 
a  very  different  relationship  with  the  elevation  rate.  During  the  translational  phase 
of  the  wing  stroke,  wing  velocity  is  only  one  component  that  influences  the  bound 
circulation  about  the  wing.  Feathering  orientation,  along  with  the  wing  velocity, 
determines  the  geometric  angle  of  attack.  Examining  the  relationship  between  the 
force  magnitude  and  the  feathering  orientation  in  Figure  33  shows  that  the  peak 
force  production  and  maximum  feathering  orientation  occur  in  phase  for  the  4-term 
kinematic  profile,  while  the  peak  force  magnitude  lags  behind  the  peak  feathering 
orientation  for  the  2-  and  3-term  kinematic  profiles. 

It  is  possible  that  the  phase  lag  between  the  peak  force  magnitude  and  the 
2-  and  3-term  sweep  and  elevation  rates  and  the  feathering  orientation  are  due  to 
phase  variations  between  the  kinematic  parameters  within  a  given  profile.  It  is  also 
possible  that  the  phase  variation  between  the  forces  and  the  parameters  are  due  to 
the  influence  of  unsteady  aerodynamic  mechanisms  that  may  not  be  represented  as 
functions  of  the  kinematic  parameters  alone.  Such  aerodynamic  mechanisms  include 
the  influence  of  interactions  between  the  wing  and  previously  shed  wake  structures, 
and  the  influence  of  local  variations  in  the  aerodynamic  mechanism  along  the  wing 
span. 
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The  4-term  kinematic  profile  exhibits  a  second  peak  in  the  force  history  as  the 
sweep  rate  is  decreasing  and  the  wing  starts  the  rotational  phase  of  motion  charac¬ 
terized  by  the  large  feathering  rotation  at  the  end  of  the  half-stroke.  The  4-term 
kinematic  profile  differers  significantly  from  the  2-  and  3-term  kinematic  profiles  in 
the  shape  of  the  feathering  parameter  and  the  magnitude  and  timing  of  the  feather¬ 
ing  rates  and  accelerations.  The  shape  of  the  2-  and  3-term  feathering  parameters  is 
nearly  sinusoidal,  while  the  4-term  feathering  parameter  exhibits  a  delayed  feather¬ 
ing  rotation,  or  a  squared  off  feathering  profile.  The  delay  corresponds  to  a  reduced 
duration  of  the  4-term  feathering  rotation,  leading  to  increased  feathering  rates  and 
accelerations.  This  variation  is  significant  because  the  second  peak  in  the  force  history 
associated  with  the  4-term  kinematic  profile  corresponds  to  the  second  peak  in  the 
feathering  acceleration  during  each  half-stroke,  as  shown  in  Figure  35.  The  influence 
of  the  feathering  acceleration  on  the  aerodynamic  mechanisms  is  examined  in  the  next 
chapter. 

Why  the  force  history  produced  by  the  4-term  kinematic  profile  exhibits  the  two 
peaks  during  each  half  stroke  is  not  yet  fully  clear.  Examining  the  phase  plots,  the 
first  aerodynamic  peak  corresponds  to  peak  translational  velocities  and  feathering  ori¬ 
entation  during  the  translational  phase  of  wing  motion,  and  the  second  aerodynamic 
peak  corresponds  to  peak  feathering  acceleration  as  the  wing  enters  the  rotational 
phase  of  motion.  Yet,  it  is  not  clear  what  causes  the  dramatic  drop  in  the  force 
magnitude  between  the  peak  translational  velocity  and  the  peak  feathering  accelera¬ 
tion.  This  phenomenon  may  also  be  observed  in  simulations  previously  computed  by 
Aono  [24,25]  for  his  examination  of  the  wake  of  a  hovering  Hawkmoth.  The  influence 
of  the  sharp  decline  in  force  magnitude  may  be  seen  as  a  dramatic  drop  at  a  constant 
feathering  rate  in  Figure  34,  and  as  dips  near  the  peak  elevation  accelerations  in  Fig¬ 
ure  35.  Understanding  this  interaction  and  identifying  a  cause  and  effect  relationship 
would  benefit  from  further  study  of  the  influence  of  the  kinematic  parameters  based 
on  a  systematic  variation  of  the  kinematic  parameters. 
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Figure  32:  Aerodynamic  force  and  moment  histories  for  the  third  flapping  cycle  of 
the  rigid  wing  solutions  based  on  the  2-,  3-,  and  4-term  Fourier  kinematic  profiles. 
The  flier  is  oriented  in  a  reference  system  with  the  cc-axis  aligned  in  the  aft  direction, 
the  £-axis  aligned  in  the  vertical  direction,  and  the  y- axis  aligned  with  the  right  wing. 
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Figure  33:  Plotting  aerodynamic  forces  against  the  time-dependent  kinematic  parameters  for  the  2-,  3-,  and  4-term  Fourier 
approximations  highlight  the  sensitivity  of  the  aerodynamic  forces  to  the  kinematic  parameters.  Rigid  wing  solution. 
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Figure  34:  Plotting  aerodynamic  forces  against  the  time-dependent  kinematic  rates  for  the  2-,  3-,  and  4-term  Fourier 

approximations  highlight  the  sensitivity  of  the  aerodynamic  forces  to  the  kinematic  rates.  Rigid  wing  solution. 
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Figure  35:  Plotting  aerodynamic  forces  against  the  time-dependent  kinematic  accelerations  for  the  2-,  3-,  and  4-term 

Fourier  approximations  highlight  the  sensitivity  of  the  aerodynamic  forces  to  the  kinematic  accelerations.  Rigid  wing 
solution. 


5.2.3  Relationship  of  the  Aerodynamic  Forces  to  Wing  Orientation.  The 
timing  of  the  maximum  and  minimum  x -,  y-,  and  ^-direction  forces  is  determined 
by  the  wing  orientation  during  the  wing  stroke  cycle.  The  peak  aerodynamic  forces 
occur  near  the  maximum  sweep  rate,  and  the  maximum  sweep  rate  occurs  at  small 
values  of  the  sweep  parameter.  At  small  sweep  angles,  the  wing  surface  normals  are 
primarily  aligned  in  the  x-z  plane,  producing  peak  x-  and  ^-direction  forces. 

Peak  side  forces  for  the  single  wing  simulations  occur  at  the  end  of  each  half¬ 
stroke  during  the  rotational  phase  of  motion.  The  rotational  phase  of  motion  is 
characterized  by  large  feathering  rotations,  during  which  the  wing  surface  normals 
are  rotating  through  the  x-y  plane.  Because  the  surface  normals  are  aligned  close  to 
the  x-y  plane,  the  pressure-based  force  contribution  in  the  side,  or  ^/-direction  is  much 
greater  than  during  the  translational  phase  of  motion. 

The  rotational  phase  of  motion  is  also  characterized  by  the  role  reversal  of  the 
“upper”  and  “lower”  wing  surfaces  and  the  associated  transition  of  the  leading  edge 
vortex  from  one  wing  surface  to  the  opposite  surface.  There  are  also  significant  non- 
circulatory  aerodynamic  forces  associated  with  the  wing  wake  interaction  and  the 
acceleration  of  the  wing  as  it  changes  direction  of  motion.  The  “upper”  and  “lower” 
wing  surface  transition  is  accompanied  by  a  sharp  change  in  the  sign  of  the  ^/-direction 
force  at  the  beginning  of  the  next  half-stroke. 

5.2.4  Spanwise  Distribution  of  Aerodynamic  Forces.  Examining  the  span- 
wise  force  distribution  highlights  spatial  variations  in  the  force  history  that  are  not 
readily  apparent  when  examining  the  time  history  of  the  net  aerodynamic  forces  and 
provides  insight  into  the  timing  and  phase  of  the  aerodynamic  forces  relative  to  the 
wing  motion.  Spatial  force  variations  also  highlight  regions  of  the  wing  influenced  by 
local  fluctuations  in  the  aerodynamic  phenomena.  Plotting  the  aerodynamic  forces 
as  a  function  of  span  and  time  provides  a  visual  correlation  between  local  wing  ori¬ 
entation  and  the  directionality  of  the  net  aerodynamic  force. 
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Spanwise  force  distributions  for  the  rigid  wing  solutions  based  on  the  4-term 
kinematic  profile  are  shown  from  the  side,  front,  and  top  perspectives  in  Figures  36, 
37,  and  38.  Following  the  spanwise  forces  through  the  flapping  cycle,  the  4-term  forces 
peak  early  in  each  half- stroke,  as  seen  at  t/r  =  1/8  and  5/8,  followed  by  a  decrease 
in  forces  in  the  middle  of  each  half-stroke,  at  t/r  =  1/2  and  1.  The  decrease  in  forces 
is  followed  by  another  aerodynamic  peak  late  in  each  half-stroke,  seen  at  t/r  =  3/8 
and  7/8.  The  aerodynamic  forces  produced  early  in  the  half-stroke  exhibit  a  spanwise 
force  distribution  that  peaks  near  the  wing  tip.  This  corresponds  to  the  observed 
distribution  of  the  strength  and  size  of  the  leading  edge  vortex,  which  produces  a 
conical  structure  increasing  in  diameter  towards  the  wing  tip.  However,  the  increase 
in  aerodynamic  forces  late  in  the  half-stroke  exhibit  a  more  uniform  force  distribution 
when  compared  to  the  force  distribution  associated  with  the  first  aerodynamic  peak. 
Compared  to  the  early  aerodynamic  peak,  the  later  peak  at  t/r  =  3/8  and  7/8  exhibits 
a  significant  force  contribution  near  the  wing  root  and  a  reduced  contribution  near 
the  wing  tip. 

As  discussed  previously,  the  second  aerodynamic  peak  corresponds  with  an  in¬ 
crease  in  the  feathering  acceleration,  or  the  rotational  acceleration  about  the  leading 
edge  of  the  wing.  Because  the  rectangular  wing  has  a  uniform  chord  length,  the  non- 
circulatory  aerodynamic  phenomena  associated  with  the  feathering  rotation  should 
have  a  consistent  influence  along  the  entire  length  of  the  wing.  However,  because 
of  three-dimensional  influences  due  to  the  finite  length  of  the  wing,  and  interactions 
with  the  other  aerodynamic  phenomena  and  the  non-uniform  wake,  the  actual  force 
distribution  is  non-uniform. 

Examining  the  force  distribution  at  the  transition  from  one  half-stroke  to  the 
next,  at  t/r  =  0  and  1/2,  the  aerodynamic  forces  are  acting  in  different  directions 
along  the  wing  length.  This  spanwise  phase  variation  in  the  directionality  of  the 
aerodynamic  forces  illustrates  how  circulatory  and  non-circulatory  forces  interact  to 
produce  the  complex  force  histories  exhibited  by  flapping  wing  fliers.  At  the  end  of 
the  half-stroke,  the  role  of  the  “upper”  and  “lower”  wing  surfaces  is  changing  as  the 
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wing  undergoes  the  large  feathering  rotation.  This  is  accompanied  by  the  release 
of  the  leading  edge  vortex  produced  during  the  preceding  the  half-stroke  into  the 
wake,  and  the  formation  of  a  new  leading  edge  vortex  on  the  opposite  wing  surface  as 
the  wing  changes  direction  of  motion.  The  leading  edge  vortex  is  strongest  towards 
the  wing  tip  where  the  translational  velocities  are  greatest,  while  the  non-circulatory 
aerodynamic  forces  associated  with  the  feathering  rotation  produce  a  more  uniform  set 
of  forces  along  the  wing  length.  It  should  be  noted  that  this  is  a  simplified  discussion 
that  neglects  the  influence  of  wing-wake  interactions  and  the  non-circulatory  forces 
induced  by  the  feathering  and  elevation  accelerations. 

The  spanwise  force  distributions  for  the  rigid  wing  solutions  based  on  the  3-term 
kinematic  profile,  shown  from  the  side,  front,  and  top  perspectives  in  Figures  39,  40, 
and  41,  exhibit  a  similar  spanwise  force  distribution  to  the  4-term  rigid  wing  solutions 
through  the  wing  stroke.  The  difference  in  spanwise  forces  early  in  the  down  stroke 
is  shown  in  Figure  42.  Atf/r  =  5/8,  the  aerodynamics  should  be  dominated  by  the 
leading  edge  vortex,  and  variations  in  the  aerodynamic  forces  will  be  due  to  variations 
in  the  strength  of  the  leading  edge  vortex. 

5.2.5  Three-Dimensional  Wake.  Figure  43  highlights  the  wake  structures 
produced  by  the  rigid  wing  at  pronation.  The  wake  structure  at  pronation  is  influenced 
by  the  deceleration  of  the  wing,  which  initializes  the  leading  edge  vortex  on  the 
“upper”  wing  surface  and  the  starting  vortex  at  the  trailing  edge  of  the  wing.  The 
wake  is  also  influenced  by  the  feathering  rotation,  which  influences  the  starting  vortex 
and  pressure  distribution  along  the  wing  surface.  Visualizing  the  wake  using  iso¬ 
surfaces  of  vorticity  in  Figures  44  and  45,  the  leading  edge  vortex  can  be  seen  behind 
the  leading  edge  on  the  “upper”  wing  surface,  expanding  in  diameter  from  the  wing 
root  to  the  wing  tip  for  the  rigid  wing  simulations  based  on  the  3-  and  4-term  kinematic 
profiles.  Other  vortical  structures,  such  at  the  tip  vortex  and  the  starting  and  stopping 
vortices,  can  also  be  identified  in  the  wake.  The  interaction  of  the  leading  edge  vortex 
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Figure  36:  Span-wise  force  distribution  for  the  rigid  wing  based  on  the  4-term 

kinematic  profile.  Left  wing  viewed  from  the  left  side.  The  three  components  of  the 
sectional  aerodynamic  force  are  shown  in  red  for  the  x-direction,  in  green  for  the 
^-direction,  and  in  bine  for  the  ^-direction.  The  total  aerodynamic  force  is  shown  in 
cyan  at  the  wing  root.  The  wing  surface  with  chordwise  marking  is  shown  in  gray 
and  the  tip  path  is  shown  in  black. 
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Figure  37:  Span-wise  force  distribution  for  the  rigid  wing  based  on  the  4-term 

kinematic  profile.  Left  wing  viewed  from  above.  The  three  components  of  the  sectional 
aerodynamic  force  are  shown  in  red  for  the  direction,  in  green  for  the  ^/-direction, 
and  in  blue  for  the  ^-direction.  The  total  aerodynamic  force  is  shown  in  cyan  at  the 
wing  root.  The  wing  surface  with  chordwise  marking  is  shown  in  gray  and  the  tip 
path  is  shown  in  black. 
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Figure  38:  Span-wise  force  distribution  for  the  rigid  wing  based  on  the  4-term 

kinematic  profile.  Left  wing  viewed  from  the  front.  The  three  components  of  the 
sectional  aerodynamic  force  are  shown  in  red  for  the  x- direct  ion,  in  green  for  the 
^/-direction,  and  in  blue  for  the  ^-direction.  The  total  aerodynamic  force  is  shown  in 
cyan  at  the  wing  root.  The  wing  surface  with  chordwise  marking  is  shown  in  gray 
and  the  tip  path  is  shown  in  black. 
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Figure  39:  Span-wise  force  distribution  for  the  rigid  wing  based  on  the  3-term 

kinematic  profile.  Left  wing  viewed  from  the  left  side.  The  three  components  of  the 
sectional  aerodynamic  force  are  shown  in  red  for  the  x-direction,  in  green  for  the 
^-direction,  and  in  bine  for  the  ^-direction.  The  total  aerodynamic  force  is  shown  in 
cyan  at  the  wing  root.  The  wing  surface  with  chordwise  marking  is  shown  in  gray 
and  the  tip  path  is  shown  in  black. 
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Figure  40:  Span-wise  force  distribution  for  the  rigid  wing  based  on  the  3-term 

kinematic  profile.  Left  wing  viewed  from  above.  The  three  components  of  the  sectional 
aerodynamic  force  are  shown  in  red  for  the  direction,  in  green  for  the  ^/-direction, 
and  in  blue  for  the  ^-direction.  The  total  aerodynamic  force  is  shown  in  cyan  at  the 
wing  root.  The  wing  surface  with  chordwise  marking  is  shown  in  gray  and  the  tip 
path  is  shown  in  black. 


82 


t/tau 


Figure  41:  Span-wise  force  distribution  for  the  rigid  wing  based  on  the  3-term 

kinematic  profile.  Left  wing  viewed  from  the  front.  The  three  components  of  the 
sectional  aerodynamic  force  are  shown  in  red  for  the  x- direct  ion,  in  green  for  the 
^/-direction,  and  in  blue  for  the  ^-direction.  The  total  aerodynamic  force  is  shown  in 
cyan  at  the  wing  root.  The  wing  surface  with  chordwise  marking  is  shown  in  gray 
and  the  tip  path  is  shown  in  black. 
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Figure  42:  Spanwise  forces  at  t/r  =  5/8  for  the  3-term  and  4-term  rigid  wing 

solutions.  The  solutions  exhibit  similar  force  distributions  with  a  significant  variation 
in  force  magnitudes. 
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Figure  43:  Annotated  wake  structures  for  the  4-term  rigid  wing  simulations  at 
t/r  =  0.6.  The  wing  is  near  the  end  of  the  feathering  rotation  at  pronation  and 
starting  the  sweep  rotation  for  the  downstroke. 

with  the  tip  vortex  and  the  starting  and  stopping  vortices  leads  to  the  complex  wake 
structures  and  influence  the  evolution  of  the  wake  structure  with  time. 

When  the  timing  and  orientation  of  the  wake  structures  align,  a  momentum  jet 
is  produced,  providing  both  lift  and  thrust  for  hover  and  forward  flight.  At  hover, 
the  momentum  jet  is  aligned  in  the  vertical  direction,  sustaining  the  weight  of  the 
Hawkmoth.  The  momentum  jet  produced  by  the  rigid  wing  3-  and  4-term  kinematic 
simulations  are  shown  in  Figures  46  and  47  using  contours  of  the  mean  velocity  field 
computed  over  the  duration  of  one  flapping  cycle. 


5.2.6  Mean  Aerodynamic  Forces  and  Aerodynamic  Power  Requirements. 
Mean  aerodynamic  forces,  moments,  and  power  requirements  were  computed  for  the 


85 


tfr  =  4/8 


3-term 


4-term 


Figure  44:  Iso-surface  of  vorticity  magnitude  for  the  3-  and  4-term  kinematic  profiles 
shown  during  the  upstroke. 
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Figure  45:  Iso-surface  of  vorticity  magnitude  for  the  3-  and  4-term  kinematic  profiles 

shown  during  the  downstroke. 
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Figure  46:  Mean  velocity  components  normalized  by  tip  velocity  for  the  3-term 

rigid  wing  simulation.  The  three  xy-planes  are  above,  immediately  below,  and  below 
the  stroke-plane.  The  coherent  vertical  momentum  jet  can  be  seen  in  the  plots  of  the 
tc-velocity  component  and  the  velocity  magnitude. 
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Figure  47:  Mean  velocity  components  normalized  by  tip  velocity  for  the  4-term 

rigid  wing  simulation.  The  three  x?/-planes  are  above,  immediately  below,  and  below 
the  stroke-plane.  The  vertical  momentum  jet  can  be  seen  in  the  plots  of  the  w- velocity 
component  and  the  velocity  magnitude. 
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2-,  3-,  and  4-term  kinematic  profiles  in  order  to  compare  the  performance.  The  mean 
aerodynamic  forces,  moments,  and  power  requirements  for  the  rigid  wing  solutions 
are  presented  in  Table  3.  Of  the  three  kinematic  profiles,  the  4-term  profile  produced 
the  largest  mean  vertical  force,  followed  by  the  2-term  and  3-term  profiles.  However, 
the  3-term  kinematic  approximation  was  the  most  efficient  of  the  three  kinematic 
approximations  in  terms  of  mean  vertical  force  compared  to  the  mean  aerodynamic 
power  required. 


Table  3:  Computed  aerodynamic  forces,  moments,  and  power  requirements  non- 

dimensionalized  by  the  dynamic  pressure  at  the  wing  tip  for  the  2-,  3-,  and  4- Term 
kinematic  profiles,  rigid  wing  solution 


2- Term 

3- Term 

4- Term 

CFx 

1.165xlO~02 

3.638x10-°° 

1.779xl0"°4 

CFy 

-2.919xl0"°2 

-1.625xl0-°2 

-2.841  xl0"°4 

CFz 

1.136xlO”01 

1.249xlO_01 

1.211xlO“01 

CMx 

-7.325xl0~°2 

-7.608  xl0-°2 

-7.390x  10-02 

CMy 

-1.388xlO~02 

-1.280xl0-°2 

-1.087x  10“°2 

CMz 

6.108x10-°° 

4.001  xl0-°3 

4.942x10-°° 

CP 

6.755xl0~°2 

7.229  xlO“02 

7.752  xlO”02 

Vaero 

1.682x10+°° 

1.728x10+°° 

1.562x10+°° 

5. 3  Summary 

Examining  the  influence  of  kinematic  variation  on  the  unsteady  aerodynamic 
phenomena  using  a  truncated  Fourier  approximation  of  the  kinematic  parameters  does 
not  allow  for  the  influence  of  a  specific  kinematic  parameter  to  be  examined.  Docu¬ 
menting  the  influence  of  a  specific  kinematic  parameter  on  aerodynamic  performance 
would  require  the  kinematic  parameters  to  be  defined  such  that  specific  parameters 
could  be  perturbed,  or  by  mixing  the  number  of  Fourier  terms  used  to  define  a  given 
kinematic  profile.  A  kinematic  profile  could  be  defined  using  a  1-term  Fourier  ap¬ 
proximation  of  elevation,  a  3-term  Fourier  approximation  of  feathering,  and  a  2-term 
Fourier  approximation  of  sweep,  and  so  forth.  However,  such  analysis  would  require 
a  significant  increase  in  the  number  of  required  simulations. 
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Examining  the  influence  of  the  Fourier  approximation  on  aerodynamic  perfor¬ 
mance  provides  insight  into  the  sensitivity  of  the  aerodynamic  forces  and  unsteady 
aerodynamic  phenomena  to  the  kinematic  approximation  used  as  a  baseline  for  exam¬ 
ining  the  influence  of  structural  response.  This  examination  highlighted  the  coupling 
between  the  strength  of  the  leading  edge  vortex,  a  circulatory  aerodynamic  mecha¬ 
nism,  with  sweep  rate  and  feathering  orientation.  Also  highlighted  was  the  influence  of 
the  feathering  acceleration  on  the  non-circulatory  aerodynamic  forces  produced  during 
the  rotational  motion  at  the  end  of  each  half-stroke.  Identifying  and  understanding 
each  of  these  interactions  provides  insight  into  how  structural  response  may  alter  the 
strength  and  timing  of  the  aerodynamic  phenomena  and  aerodynamic  forces,  as  the 
influence  of  structural  response  may  be  conceptualized  as  local  kinematic  variations. 

Willmott  and  Ellington  noted  that  the  4-term  Fourier  approximation  of  the  kine¬ 
matic  parameters  defined  in  Equations  (22),  (23),  and  (24)  adequately  represented  the 
experimentally  documented  wing  motion  while  minimizing  the  influence  of  spurious 
higher  order  harmonics  due  to  the  experimental  approach  on  the  kinematic  rate  and 
acceleration  profiles  [6].  Comparing  the  aerodynamic  force  histories  and  the  unsteady 
aerodynamic  phenomena  produced  by  the  2-,  3-,  and  4-term  kinematic  profiles  pre¬ 
sented  in  this  chapter,  the  aerodynamic  response  is  clearly  sensitive  to  variations  in 
the  kinematic  rate  and  acceleration  profiles,  as  well  as  phase  variations  between  the 
kinematic  parameters. 
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VI.  Analysis  of  Feathering  Variation 

As  discussed  in  the  previous  chapter,  the  unsteady  aerodynamic  phenomena  are  sen¬ 
sitive  to  variations  in  velocity,  acceleration,  and  orientation  of  the  wing  throughout 
the  wing  stroke.  It  has  been  demonstrated  that  variations  in  the  rigid  wing  kine¬ 
matics  impact  the  time  history  of  the  aerodynamic  forces  and  the  wake  structures 
by  altering  the  local  aerodynamic  phenomena.  Incorporating  structural  response  into 
the  current  analysis  introduces  local  variations  into  the  velocity  profile  and  orienta¬ 
tion  of  the  wing.  These  local  variations  have  the  potential  to  influence  aerodynamic 
performance  in  a  similar  manner  to  the  kinematic  variations.  Yet,  the  influence  of 
structural  response  on  aerodynamic  performance  has  not  been  well  studied. 

To  understand  the  influence  of  structural  flexibility  on  the  underlying  unsteady 
aerodynamic  mechanisms  and  aerodynamic  performance,  deformations  of  the  wing 
surface  were  prescribed  as  a  function  of  spanwise  location  and  time  in  the  wing  local 
reference  frame.  When  these  local  deformations  are  superimposed  onto  the  rigid  wing 
motion  defined  by  the  kinematic  parameters,  it  is  possible  to  mimic  the  flexible  wing 
motion  observed  for  natural  fliers  in  free-flight. 

OVERFLOW  2.1-Elastic  is  able  to  simulate  the  flow-held  about  flexible  bodies 
exhibiting  “beam  type”  deformations  where  the  time-dependent  deformation  of  the 
body  is  described  using  a  combination  of  local  displacements,  Sx,  Sy,  Sz,  and  rota¬ 
tions,  9X,  9y,  9Z,  along  a  defined  elastic  axis.  Applied  to  the  happing  wing  problem, 
these  beam-type  deformations  may  be  used  to  model  feathering,  or  torsion  about  the 
feathering  axis,  and  bending  along  the  wing  length.  The  algebraic  mesh  deformation 
scheme  is  not  able  to  accommodate  “plate  type”  deformations  associated  with  chord- 
wise  deformations,  or  adaptive  camber,  of  the  wing  section  as  implemented.  The 
types  of  hexible  wing  motion  that  may  be  modeled  using  this  approach  are  illustrated 
in  Figure  48. 

With  the  ability  to  simulate  hexible  wing  structures  using  local  rotations  and 
deformations  along  a  dehned  elastic  axis,  the  limited  kinematic  data  documenting  the 
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spanwise  feathering  variation  was  used  to  examine  the  influence  of  spanwise  feathering 
on  the  local  aerodynamic  phenomena. 

6.1  Prescribed  Feathering  as  a  Function  of  Span  and  Time 

The  feathering  profile  used  by  Liu  [19]  and  Aono  [24]  to  define  the  “rigid-wing” 
motion  of  the  Hawkmoth  wing  in  time  is  based  on  a  mean  of  the  feathering  profiles 
recorded  at  multiple  spanwise  locations.  The  rigid  feathering  profile  is  based  on 
feathering  profiles  at  spanwise  locations  from  r/R  —  0.6  to  0.9  recorded  by  Willmott 
and  Ellington  [17].  The  feathering  profiles  measured  along  the  wing  length  form  the 
basis  of  the  simulations  incorporating  prescribed  feathering  variations. 

6.1.1  The  Flexible  Feathering  Parameter.  As  with  the  time  history  of  the 
kinematic  parameters,  the  spatial  distribution  of  the  feathering  parameter  must  be 
smoothed  in  such  a  way  that  the  feathering  variation  is  continuous  from  the  wing  root 
to  wing  tip  for  use  with  the  numerical  simulations.  The  spatial  distribution  of  the 
feathering  parameter  could  be  represented  using  a  set  of  experimentally  determined 
modal  contributions  from  a  characterization  of  the  Hawkmoth  wing  structure.  Lack¬ 
ing  this  structural  characterization  of  the  Hawkmoth  wing,  it  was  assumed  for  this 
research  that  the  spanwise  feathering  distribution  could  be  represented  using  the  first 
torsion  mode  of  a  fixed-free  beam  defined  in  Equation  (25). 


The  full  definition  of  the  “flexible”  feathering  parameter  in  space  and  time, 
a(r,  t),  depends  on  a  spatial  distribution  of  the  feathering  parameter  and  time-dependent 
boundary  conditions  to  scale  the  feathering  variation  in  time.  The  amplitude  of  the 
first  torsion  mode  defined  in  Equation  (25)  ranges  from  0  to  1  along  the  wing  length. 

If  the  orientation  of  the  wing  root  is  fixed,  or  constant  in  time,  and  the  feathering 
orientation  at  the  wing  tip  is  centered  about  zero,  the  spatial  distribution  may  simply 
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Figure  48:  Examples  of  the  flexible  wing  motion  during  the  upstroke  incorporating 

bending  and  torsion  modes  for  an  elastic  axis  along  the  leading  edge  of  the  wing.  The 
red  line  represents  the  flexible  wing  tip  path  during  the  wing  stroke  and  the  blue  line 
represents  the  rigid  wing  tip  path.  These  are  the  types  of  deformations  that  may  be 
captured  using  OVERFLOW  2.1-Elastic. 
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be  scaled  by  the  amplitude  of  the  feathering  orientation  at  the  wing  tip  to  define  the 
flexible  feathering  orientation  as  a  function  of  time  and  wing  length,  as  in  Equation 
(26). 

a(r,  t)  =  a„p(t)rT(r)  (26) 

For  cases  where  the  wing  root  is  not  fixed  in  time  or  the  feathering  parameter 
is  not  centered  around  zero,  as  is  the  case  for  the  Hawkmoth  kinematics,  a  flexible 
feathering  distribution  may  be  defined  using  the  amplitude  of  the  feathering  variation 
between  the  wing  root  and  wing  tip,  5a(t)  =  atiP(t )  —  aroot(t),  and  the  feathering 
orientation  at  the  wing  root,  as  defined  in  Equation  (27). 

a(r,  t)  =  aroot(t )  +  Sa(t)VT(r)  (27) 

To  maintain  consistency  with  the  rigid  wing  simulations,  the  rigid  wing  feathering 
parameter  was  used  to  define  the  wing  tip  orientation,  while  the  amplitude  of  the 
feathering  orientation  at  the  wing  root  was  varied  to  represent  wings  with  varying 
amounts  of  torsional  flexibility. 

The  ability  of  the  first  torsion  mode  to  capture  the  spanwise  feathering  distri¬ 
bution  measured  by  Willmott  and  Ellington  [17]  is  shown  in  Figure  49.  The  first 
torsion  mode  captures  the  dominant  feathering  variation  from  the  wing  root  to  the 
wing  tip.  A  better  capture  of  the  recorded  feathering  variation  may  be  achieved  with 
the  inclusion  of  higher-order  torsion  modes  once  a  structural  characterization  of  the 
Hawkmoth  wing  is  available. 

6.2  Influence  of  the  Feathering  Variation  on  the  Aerodynamic  Response 

To  assess  the  influence  of  the  spanwise  feathering  variation  on  the  unsteady 
aerodynamic  phenomena  and  aerodynamic  performance,  simulations  based  on  the  3- 
and  4-term  kinematic  profiles  were  computed  for  three  amplitudes  of  the  feather¬ 
ing  variation,  Sa.  The  rigid  wing  simulations,  discussed  in  Chapter  V,  represent  0% 
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Figure  49:  The  flexible  feathering  parameter,  a(r,t)  approximated  using  a  first- 

torsion  span- wise  distribution  where  aroot(t)  =  90  degrees  and  aup(t )  is  based  on  the 
rigid  wing  feathering  orientation.  The  measured  feathering  angles  were  digitized  from 
Willmott  and  Ellington  [17]. 


96 


feathering  where  aroot(t )  =  atip(t).  Solutions  were  computed  for  a  wing  with  a  re¬ 
duced  amplitude  feathering  profile  at  the  wing  root,  aroot(t)  =  0.5 atiP(t),  representing 
50%  feathering,  and  for  a  wing  with  a  fixed  feathering  orientation  at  the  wing  root, 
aroot(t )  =  90  degrees,  representing  100%  feathering.  This  0%,  50%,  and  100%  feather¬ 
ing  nomenclature  is  used  throughout  the  remainder  of  this  document  when  discussing 
the  influence  of  prescribed  feathering. 

As  discussed  in  the  previous  chapter,  the  4-term  kinematic  profile  was  consid¬ 
ered  by  Willmott  and  Ellington  to  be  a  reasonable  reconstruction  of  the  Hawkmoth 
wing  motions  [6].  Though  the  3-term  and  4-term  kinematic  profiles  produce  simi¬ 
lar  tip  paths,  the  kinematic  rates  and  velocities  vary  significantly.  This  variation  in 
the  velocity  and  acceleration  profile  has  a  significant  influence  on  the  force  histories 
produced  by  the  3-  and  4-term  kinematic  profiles,  with  the  3-term  kinematic  profile 
producing  a  single  peak  in  the  vertical  force  per  half  stroke  and  the  4-term  kinematic 
profile  producing  two  peaks  in  the  vertical  force  per  half-stroke.  The  3-term  kine¬ 
matic  profile  represents  a  simplification  in  the  underlying  wing  kinematics  compared 
to  the  4-term  profile,  as  might  be  expected  from  simplified  kinematic  actuators  that 
may  used  in  the  next  generation  of  micro-aerial  vehicles.  The  use  of  both  3-  and  4- 
terrn  kinematic  profiles  for  examining  the  influence  of  structural  flexibility  allows  for 
trends  in  the  aerodynamic  response  to  be  identified  that  are  not  limited  to  a  specific 
kinematic  profile. 

6.2.1  Influence  on  the  Aerodynamic  Force  History.  Increasing  the  torsional 
flexibility  by  decreasing  the  amplitude  of  the  feathering  profile  at  the  wing  root  while 
maintaining  a  consistent  feathering  profile  at  the  wing  tip  influences  the  timing  and 
amplitude  of  the  net  aerodynamic  forces  produced  by  the  3-term  kinematic  profile, 
shown  in  Figure  50.  The  duration  and  magnitude  of  the  single  peak  in  the  force 
magnitude  and  z-force  component  produced  by  the  rigid  wing,  with  0%  feathering, 
changes  as  the  amount  of  feathering  “flexibility”  increases.  Both  the  upstroke  and 
downstroke  see  an  increase  in  the  force  magnitude  early  in  the  half-stroke,  leading 
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to  a  longer  duration  of  the  sustained  force  magnitude.  However,  the  peak  value  of 
the  2-force  component  and  force  magnitude  during  the  downstroke  decrease  with  an 
increase  in  the  feathering  variation. 

The  timing  of  the  changes  in  the  force  history  provide  clues  to  which  aero¬ 
dynamic  phenomena  are  sensitive  to  changes  in  the  spanwise  feathering  variation. 
Changes  in  the  aerodynamic  forces  during  the  early  and  middle  portions  of  the  half¬ 
stroke  influence  the  translational  phase  of  motion  that  is  dominated  by  the  forma¬ 
tion  of  the  leading  edge  vortex  and  associated  circulatory  aerodynamic  phenomena. 
Changes  in  the  aerodynamic  forces  during  pronation  and  supination  indicate  changes 
in  the  non-circulatory  aerodynamic  phenomena  influenced  by  the  feathering  acceler¬ 
ation  and  wing-wake  interactions  that  occur  as  the  wing  changes  direction  of  motion, 
as  discussed  in  Section  5.2.2. 

At  first  glance,  “increased  torsional  flexibility”  appears  to  have  a  different  influ¬ 
ence  on  the  aerodynamic  forces  produced  by  the  4-term  kinematic  profile,  shown  in 
Figure  51,  than  observed  for  the  3-term  kinematic  profile,  shown  in  Figure  50.  How¬ 
ever,  reducing  the  amplitude  of  the  root  feathering  profile  increases  the  magnitude 
of  the  initial  aerodynamic  peak  during  the  early  portion  of  each  half-stroke,  while 
decreasing  the  amplitude  and  altering  the  phase  of  the  second  aerodynamic  peak  late 
in  the  half-stroke.  This  is  a  similar  behavior  to  that  observed  for  the  3-term  kinematic 
profile  with  an  increase  in  torsional  flexibility,  if  the  underlying  differences  between 
the  3-  and  4-term  force  histories  are  neglected. 

6.2.2  Spanwise  Velocity  and  Force  Distributions.  Introducing  a  spanwise 
feathering  variation  influences  the  geometric  configuration  of  the  wing  during  both 
phases  of  wing  motion,  translation  and  rotation.  During  the  rotational  phase  of  mo¬ 
tion,  the  reduced  amplitude  of  root  feathering  orientation  reduces  the  accompanying 
feathering  rate  and  acceleration  of  the  inboard  portion  of  the  wing.  This  change  in 
the  feathering  rate  influences  the  velocity  distribution  along  the  wing  length  at  every 
point  not  located  on  the  feathering  axis.  The  translational  velocity  of  the  wing  at  the 
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Figure  50:  Aerodynamic  force  and  moment  histories  for  the  torsionally  flexible 

simulations.  The  wing  kinematics  are  based  on  the  3-term  Fourier  approximation. 
The  flier  is  oriented  in  a  reference  system  with  the  x-axis  aligned  in  the  aft  direction, 
the  £-axis  aligned  in  the  vertical  direction,  and  the  y- axis  aligned  with  the  right  wing. 
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Figure  51:  Aerodynamic  force  and  moment  histories  for  the  torsionally  flexible 

simulations.  The  wing  kinematics  are  based  on  the  4-term  Fourier  approximation. 
The  flier  is  oriented  in  a  reference  system  with  the  x-axis  aligned  in  the  aft  direction, 
the  £-axis  aligned  in  the  vertical  direction,  and  the  y- axis  aligned  with  the  right  wing. 
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three-quarter  chord  location  is  shown  in  Figures  52,  and  53  using  a  set  of  red  velocity 
vectors  for  the  rigid  wing,  0%  feathering,  and  the  100%  feathering  wing  as  the  wing 
approaches  pronation,  at  t/r  =  3/8,  and  just  after  stroke  reversal,  at  t/r  =  4/8.  Fig¬ 
ure  52  represents  the  force  and  velocity  distribution  for  the  3-term  kinematic  profile 
and  Figure  53  represents  the  force  and  velocity  distribution  for  the  4-term  kinematic 
profile. 

At  t/r  =  3/8,  the  rigid  wing  velocity  distribution  for  both  kinematic  profiles  are 
linear  at  the  three-quarter  chord  location  and  exhibits  a  non-zero  velocity  component 
at  the  wing  root.  This  velocity  distribution  is  accompanied  by  a  spanwise  force 
distribution  at  the  leading  edge  of  the  wing,  shown  using  a  set  of  blue  vectors  in 
Figures  52  and  53,  that  exhibits  a  significant  force  contribution  near  the  wing  root. 
Contrasting  the  rigid  wing  velocity  distribution  to  the  torsionally  flexible  wing  with 
100%  feathering  variation,  the  velocity  distribution  of  the  100%  feathering  wing  at  the 
three-quarter  chord  is  not  linear,  but  mirrors  the  spatial  distribution  of  the  feathering 
parameter.  Since  the  orientation  of  the  wing  root  is  fixed  at  aroot(t )  =  90  degrees 
for  the  100%  feathering  wing,  the  translational  velocity  of  the  wing  approaches  zero 
at  the  wing  root.  This  change  in  the  velocity  distribution  from  the  rigid  wing  is 
accompanied  by  a  change  in  the  spanwise  force  distribution.  The  force  distribution 
of  the  100%  feathering  wing  exhibits  a  decreased  contribution  near  the  wing  root, 
which  corresponds  to  the  decrease  in  the  velocity  distribution  along  the  three  quarter 
chord  near  the  wing  root.  The  middle  third  of  the  wing  exhibits  an  increase  in  force 
amplitude  when  compared  to  the  0%  feathering  wing. 

The  influence  of  the  spanwise  variation  in  feathering  rates  is  even  more  pro¬ 
nounced  when  contrasting  the  force  distribution  between  the  rigid  and  feathering 
wings  just  after  stroke  reversal,  at  t/r  =  4/8,  also  shown  in  Figure  52  for  the  3-term 
kinematic  profile  and  Figure  53  for  the  4-term  kinematic  profile.  At  t/r  =  4/8,  the 
wing  has  already  changed  direction  of  motion,  as  shown  by  the  single  red  vector  at 
the  wing  tip  along  the  leading  edge  of  the  wing.  This  location  is  coincident  with 
the  feathering  axis  and  not  influenced  by  the  feathering  rate.  However,  because  the 
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wing  is  still  pronating  at  the  start  of  the  downstroke,  the  contribution  of  the  sweep 
rate  to  the  translational  velocity  at  the  wing  tip  is  negated  by  the  contribution  of 
the  feathering  rate  about  the  leading  edge  at  the  three  quarter  chord.  For  the  rigid 
wing,  this  leads  to  an  “inverse”  velocity  distribution  along  the  three-quarter  chord 
of  the  wing  with  a  greater  translational  velocity  at  the  wing  root  than  the  wing  tip. 
This  “inverse”  velocity  distribution  along  the  three-quarter  chord  is  accompanied  by 
a  spanwise  force  distribution  that  shows  significant  forces  on  the  “lower”  wing  surface 
near  the  wing  root  and  the  “upper”  wing  surface  near  the  wing  tip.  Compared  to  the 
0%  feathering  wing,  the  100%  feathering  wing  shows  a  velocity  distribution  along  the 
three-quarter  chord  that  is  nearly  zero  along  the  length  of  the  wing.  This  is  accom¬ 
panied  by  a  spanwise  force  distribution  that  does  not  exhibit  “reversed”  forces  along 
the  length  of  the  wing. 

The  influence  of  the  altered  force  distributions  on  the  net  aerodynamic  forces 
may  be  seen  in  the  force  histories  at  t/r  =  3/8  and  t/r  =  4/8,  shown  in  Figures  50 
and  51  for  the  3-  and  4-term  kinematic  profiles.  Due  to  the  orientation  of  the  wing 
surface,  the  change  in  the  force  distribution  primarily  influences  the  net  aerodynamic 
forces  in  the  x-  and  //-directions.  The  100%  feathering  wing  shows  an  increase  in  net 
x-  and  //-forces  compared  to  the  0%  feathering  wing,  or  rigid  wing,  at  both  t/r  =  3/8 
and  t/r  =  4/8.  The  influence  is  most  significant  at  t/r  =  4/8  because  the  “reversed” 
forces  of  the  rigid  wing  act  in  opposing  directions,  reducing  the  net  aerodynamic 
forces. 

During  the  translational  phase  of  motion,  the  feathering  variation  influences  the 
orientation  of  the  wing  relative  to  the  direction  of  motion  at  the  leading  edge  due  to 
the  sweep  and  elevation  rates.  This  change  in  wing  orientation,  along  with  the  change 
in  wake  structure  near  the  wing  root,  to  be  discussed  in  the  next  section,  influences 
the  spanwise  force  distribution  by  increasing  the  amplitude  of  the  aerodynamic  forces 
in  the  inner  and  middle  thirds  of  the  wing  length  at  t/r  =  1/8  and  t/r  =  2/8,  shown 
in  Figure  54  for  the  3-term  kinematic  profile  and  Figure  55  for  the  4-term  kinematic 
profile. 
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Figure  52:  Comparison  of  the  spanwise  force  distribution  in  blue  to  the  linear 

velocity  of  the  wing  at  the  three-quater  chord  location  in  red  for  the  rigid  wing  and 
prescribed  feathering  solutions.  The  wing  motion  is  based  on  the  3-term  kinematic 
profile. 
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Figure  53:  Comparison  of  the  spanwise  force  distribution  in  blue  to  the  linear 

velocity  of  the  wing  at  the  three-quater  chord  location  in  red  for  the  rigid  wing  and 
prescribed  feathering  solutions.  The  wing  motion  is  based  on  the  4-term  kinematic 
profile. 
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Figure  54:  Comparison  of  the  spanwise  force  distribution  in  blue  to  the  linear 

velocity  of  the  wing  at  the  three-quater  chord  location  in  red  for  the  rigid  wing  and 
prescribed  feathering  solutions.  The  wing  motion  is  based  on  the  3-term  kinematic 
profile. 
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Figure  55:  Comparison  of  the  spanwise  force  distribution  in  blue  to  the  linear 

velocity  of  the  wing  at  the  three-quater  chord  location  in  red  for  the  rigid  wing  and 
prescribed  feathering  solutions.  The  wing  motion  is  based  on  the  4-term  kinematic 
profile. 
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The  point  of  view  at  t/r  =  1/8  and  t/r  =  2/8  mask  the  spanwise  velocity 
distribution  below  the  wing  surface  for  each  case  except  for  the  rigid  wing  at  t/r  =  1/8. 
The  spanwise  velocity  distribution  at  the  three  quarter  chord  exhibits  a  “reversed” 
distribution  where  the  three  quarter  chord  near  the  wing  root  is  moving  to  the  left, 
or  forward,  while  the  wing  tip  is  translating  to  the  right,  or  aft.  As  with  the  rigid 
wing  during  pronation  at  t/r  =  4/8,  this  “reversed”  velocity  distribution  exhibits  a 
decreased  force  contribution  near  the  wing  root  when  compared  to  the  100%  feathering 
simulation. 

6.2.3  Pressure  Distribution  and  Wake  Structures.  The  change  in  spanwise 
force  distribution  with  increased  feathering  variation  has  a  direct  connection  to  the 
pressure  distribution  along  the  wing  surface  and  the  wake  structures  produced  by  the 
unsteady  aerodynamic  phenomena.  One  example  of  this  connection  may  be  seen  by 
comparing  the  spanwise  force  distribution  at  pronation,  t/r  —  4/8,  in  Figure  52  to  the 
vortical  wake  structures  and  pressure  distribution  along  the  wing  surface  at  pronation 
shown  in  Figure  56  for  the  3-term  kinematic  profile. 

Contrasting  the  spanwise  force  distribution  of  the  0%  feathering  wing  with  the 
100%  feathering  wing,  the  inboard  portion  of  the  0%  feathering  wing  produces  signif¬ 
icantly  smaller  force  contribution  than  the  100%  feathering  wing  during  the  transla¬ 
tional  phase  of  motion,  and  a  “reversed”  force  contribution  during  the  peak  feathering 
rotation  at  the  end  of  the  half-stroke.  Examining  Figure  56,  the  two  feathering  profiles 
produce  different  types  of  wake  structures  near  the  wing  root  during  the  feathering 
rotation  at  t/r  =  4/8.  The  vortical  structure  produced  by  the  0%  feathering  wing, 
or  rigid  wing,  wraps  around  the  wing  root  to  influence  both  the  “upper”  and  “lower” 
wing  surface.  The  100%  feathering  wing  profile  does  not  produce  a  significant  vortical 
structure  at  the  wing  root,  rather  there  is  a  vortical  structure  that  starts  at  the  trailing 
edge  of  the  wing  and  extends  into  the  flowheld  during  the  half-stroke.  This  vortical 
structure  shed  by  the  100%  feathering  wing  “connects”  with  the  vortical  structure 
shed  from  the  wing  tip,  forming  a  ring  type  structure  at  the  end  of  the  half-stroke. 
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Figure  56:  Iso-surface  of  vorticity  magnitude  for  the  0%  and  100%  Feathering 

simulations  based  on  the  3-term  kinematic  profiles  shown  during  pronation. 
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This  structure  then  interacts  with  the  “starting”  and  “stopping”  vortical  structures 
shed  from  the  trailing  edge  of  the  wing  as  the  wing  changes  direction  of  motion  at  the 
end  of  each  half-stroke  [16].  Components  of  these  structures  may  be  identified  in  the 
wake  of  the  0%  feathering  profile.  However,  the  vortical  structures  produced  by  the 
rigid  wing  do  not  exhibit  the  same  continuous  vortical  structures  as  those  produced 
by  the  100%  feathering  wing. 

The  influence  of  the  altered  wake  structures  on  the  pressure  distribution  may  be 
seen  in  the  mean  pressure  distribution  computed  during  the  upstroke  and  downstroke 
shown  in  Figures  57  and  58  for  the  0%  and  100%  feathering  profiles.  Each  figure  shows 
two  wing  pairs.  The  wing  pair  on  the  right  side  of  the  figure  represents  the  mean 
pressure  coefficients  computed  over  the  duration  of  the  upstroke,  and  the  wing  pair 
on  the  left  side  represents  the  mean  pressure  coefficient  computed  over  the  duration 
of  the  downstroke.  The  wing  pairs  are  oriented  as  they  would  be  at  the  end  of 
their  respective  half-strokes.  Because  the  end  of  each  half-stroke  occurs  during  the 
feathering  rotation,  the  point  of  view  shows  the  “upper”  wing  surface  in  the  near-field 
and  the  “lower”  wing  surface  in  the  far-held. 

Contrasting  the  mean  pressure  coefficients  for  the  0%  and  100%  feathering  sim¬ 
ulations,  the  rigid  wings  produce  a  significant  region  of  low  pressure,  or  negative 
pressure  coefficient,  on  the  “lower”  wing  surface  during  the  upstroke,  as  seen  by  the 
green  region  on  the  “lower”  wing  surface  in  Figure  57.  This  is  accompanied  by  a 
similar,  but  smaller,  region  of  low  pressure  on  the  “lower”  wing  surface  during  the 
downstroke.  Neither  the  upstroke  or  the  downstroke  of  the  100%  feathering  case 
exhibit  these  regions  of  low  pressure  on  the  lower  wing  surface,  shown  in  Figure  58. 

The  larger  region  of  reduced  pressure  on  the  lower  wing  surface  during  the 
upstroke  corresponds  to  the  longer  delay  between  the  start  of  the  upstroke  and  the 
peak  aerodynamic  force  magnitude  than  observed  during  the  downstroke  for  the  rigid 
wing  simulation  in  Figure  50. 
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Figure  57:  Mean  preasure  coefficients  on  the  wing  surface,  3-term  kinematic  profile, 
rigid  wing  solution. 

Along  with  the  influence  of  the  feathering  variation  on  the  “lower”  wing  surface 
pressures  during  each  half-stroke,  there  is  a  corresponding  smaller,  but  still  significant, 
influence  on  the  “upper”  wing  surface  pressures.  The  mean  pressure  distribution 
along  the  “upper”  wing  surface  varies  near  the  wing  root  and  near  the  trailing  edge 
of  the  wing  between  the  0%  and  100%  feathering  simulations.  These  regions  are 
influenced  by  the  change  in  feathering  orientation,  and  by  the  interaction  of  the 
unsteady  aerodynamic  phenomena. 

A  similar  correlation  between  the  wake  structures,  Figure  59,  and  the  mean 
pressure  coefficients,  shown  in  Figures  60  and  61,  holds  for  the  4-term  kinematic 
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Figure  58:  Mean  preasure  coefficients  on  the  wing  surface,  3-term  kinematic  profile, 
100%  feathering  solution. 
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profile  as  was  described  for  the  3-term  kinematic  profile.  The  rigid  wing  exhibits 
regions  of  negative  pressure  coefficient  on  the  “lower”  wing  surface  near  the  wing  root, 
with  a  bias  in  the  influence  on  the  upstroke  compared  to  the  downstroke.  This  matches 
with  the  variation  in  the  vortical  wake  structures  observed  between  the  0%  and  100% 
feathering  cases  for  the  4-term  kinematic  profile  in  Figures  62  and  63. 

6.2.4  Wake  Velocities  and  the  Momentum  Jet.  The  vortical  wake  structures 
produced  by  the  dynamic  wing  motions  are  continually  evolving  in  both  space  and 
time.  Dominant  flow  features,  such  as  the  leading  edge  vortex  and  tip  vortex,  may  be 
identified  and  tracked  in  order  to  compare  the  wake  structures  produced  by  different 
wing  kinematics  of  feathering  profiles,  but  these  dominant  wake  structures  do  not 
provide  much  insight  into  the  strength  or  directionality  of  the  induced  momentum  jet 
produced  by  the  dynamic  wing  motion. 

Examining  the  mean  velocity  held  computed  over  a  complete  happing  cycle 
highlights  the  strength,  directionality,  and  “coherence”  of  the  vertical  momentum 
jet  produced  by  the  hover  kinematics.  Examining  the  velocity  held  for  the  4-term 
kinematic  prohle  shown  in  Figures  64  and  65,  the  wake  structures  produced  near  the 
wing  root  by  the  rigid  wing  introduce  perturbations  in  u  and  v  that  decrease  the 
“coherence”  of  the  wake  structures  as  they  evolve  in  time  compared  to  the  100% 
feathering  solution.  These  perturbations  influence  the  strength  and  direction  of  the 
induced  momentum  jet,  which  influences  aerodynamic  performance. 

6.3  Change  in  the  Phase  Relationship  Between  the  Aerodynamic  Forces 

and  Kinematic  Parameters 

The  change  in  force  history  with  increasing  feathering  variation  alters  the  phase 
relationship  between  the  aerodynamic  forces  and  the  underlying  kinematic  param¬ 
eters,  rates,  and  accelerations.  As  discussed  in  Chapter  V,  the  aerodynamic  forces 
for  the  3-  and  4-term  kinematic  profiles  are  sensitive  to  the  sweep  and  elevation  rates 
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Figure  59:  Iso-surface  of  vorticity  magnitude  for  the  0%  and  100%  Feathering 

simulations  based  on  the  4-term  kinematic  profiles  shown  during  pronation. 
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Figure  60:  Mean  preasure  coefficients  on  the  wing  surface,  4-term  kinematic  profile, 
rigid  wing  solution. 
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Figure  61:  Mean  preasure  coefficients  on  the  wing  surface,  4-term  kinematic  profile, 
100%  feathering  solution. 
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Figure  62:  Iso-surface  of  vorticity  magnitude  for  the  0%  and  100%  Feathering 

simulations  based  on  the  4-term  kinematic  profiles  shown  during  the  upstroke. 
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Figure  63:  Iso-surface  of  vorticity  magnitude  for  the  0%  and  100%  Feathering 

simulations  based  on  the  4-term  kinematic  profiles  shown  during  the  downstroke. 


117 


-3 


-3 


X 10  X 10 


Figure  64:  Mean  velocity  components  normalized  by  tip  velocity  for  the  4-term 

rigid  wing  simulation.  The  three  x?/-planes  are  above,  immediately  below,  and  below 
the  stroke-plane.  The  vertical  momentum  jet  can  be  seen  in  the  plots  of  the  w- velocity 
component  and  the  velocity  magnitude. 
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Figure  65:  Mean  velocity  components  normalized  by  tip  velocity  for  the  4-term 

rigid  wing  simulation.  The  three  x?/-planes  are  above,  immediately  below,  and  below 
the  stroke-plane.  The  vertical  momentum  jet  can  be  seen  in  the  plots  of  the  w- velocity 
component  and  the  velocity  magnitude. 
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and  feathering  orientation  during  the  translational  phase  of  motion,  and  the  feathering 
acceleration  during  the  rotational  phase  of  motion. 

Examining  the  phase  relationship  for  the  4-term  kinematic  profile  shown  in  Fig¬ 
ures  66,  67,  and  68,  it  appears  that  increasing  the  feathering  variation  by  decreasing 
the  amplitude  of  the  root  feathering  profile  heightens  the  “correlation”  between  the 
aerodynamic  force  and  some  kinematic  parameters  while  reducing  the  influence  of 
other  parameters.  However,  incorporating  the  feathering  variation  does  not  funda¬ 
mentally  change  the  phase  relationship  exhibited  by  any  of  the  kinematic  parameters. 

Increasing  the  feathering  variation  heightens  the  phase  relationship  to  sweep 
rate,  elevation  rate,  and  feathering  orientation  during  the  translational  phase  of  mo¬ 
tion,  as  can  be  seen  by  the  sharp  peaks  in  the  phase  portraits  exhibited  by  the  100% 
feathering  case.  This  is  accompanied  by  an  increase  in  the  sensitivity  to  feathering 
acceleration  rate  for  the  100%  feathering  case  during  pronation  and  supination.  One 
surprise  influence  is  the  heightened  sensitivity  of  the  aerodynamic  forces  to  the  eleva¬ 
tion  parameter,  shown  in  Figure  66,  and  elevation  acceleration,  shown  in  Figure  68. 
The  increase  in  duration  of  the  aerodynamic  force  magnitude  caused  by  the  earlier 
rise  in  forces  during  the  half-stroke  alters  the  phase  relationship  between  the  sweep 
parameter,  elevation  rate,  feathering  rate,  and  sweep  acceleration  as  the  feathering 
variation  increases. 

6.4  Influence  on  Aerodynamic  Performance 

Comparing  the  simulations  on  the  basis  of  aerodynamic  efficiency,  the  perfor¬ 
mance  metrics  show  that  the  3-term  kinematic  profile,  Table  4,  was  the  most  efficient, 
followed  closely  by  the  4-term  kinematic  profile,  Table  5.  Comparing  the  simulations 
on  the  basis  of  aerodynamic  power  requirements,  the  3-term  rigid  wing  simulation  re¬ 
quired  the  least  amount  of  aerodynamic  power,  requiring  37%  less  aerodynamic  power 
than  the  4-term  100%  feathering  simulation. 
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Figure  66:  Aerodynamic  force  magnitude  as  a  function  of  kinematic  parameter  for  a  rigid  wing,  0%,  half  amplitude  root 

feathering  profile,  50%,  and  a  fixed  root  feathering  orientation,  100%,  based  on  the  4-term  kinematic  profile. 
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Figure  67:  Aerodynamic  force  magnitude  as  a  function  of  kinematic  rate  for  a  rigid  wing,  0%,  half  amplitude  root 

feathering  profile,  50%,  and  a  fixed  root  feathering  orientation,  100%,  based  on  the  4-term  kinematic  profile. 
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Figure  68:  Aerodynamic  force  magnitude  as  a  function  of  kinematic  acceleration  for  a  rigid  wing,  0%,  half  amplitude 

root  feathering  profile,  50%,  and  a  fixed  root  feathering  orientation,  100%,  based  on  the  4-term  kinematic  profile. 


Within  the  3-term  and  4-term  kinematic  profiles,  the  solution  computed  for  a 
root  feathering  amplitude  of  50%  of  the  wing  tip  feathering  amplitude  generated  the 
most  vertical  force  and  was  the  most  efficient  within  a  given  kinematic  profile. 

In  addition  to  peak  vertical  force  and  power  requirements,  the  solutions  with  pre¬ 
scribed  feathering  variations  produce  significantly  smaller  mean  aerodynamic  forces 
in  the  side  and  fore-aft  directions,  the  y-  and  x- direct  ions  respectively.  The  reduction 
in  mean  fore-aft  forces  ranged  from  a  factor  of  2  to  a  factor  of  65,  and  the  reduction 
in  mean  side  forces  ranged  from  a  factor  of  2  to  a  factor  of  100,  which  has  significant 
implications  for  future  work  in  guidance  and  control  of  engineered  flapping  wing  fliers. 


Table  4:  Computed  aerodynamic  forces,  moments,  and  power  for  prescribed  feath¬ 
ering  variation  using  the  3-term  kinematic  profile 


0%  Feathering 

50%  Feathering 

100%  Feathering 

CFx 

1.164xlO~02 

3.638  xlO”03 

1.779  xlO"04 

CFy 

-2.919xlO”02 

-1.624xlO“02 

-2.840  xlO"04 

CFz 

1.136xlO”01 

1.249xlO“01 

1.210xl0“01 

CMx 

-7.325xlO“02 

-7.607xl0-02 

-7.390  xlO~02 

CMy 

-1.387xlO~02 

-1.280xl0-02 

-1.086xl0-02 

CMz 

6.108xl0"°3 

4.000  xlO”03 

4.941  xlO’03 

CP 

6.755xlO-02 

7.229  xlO”02 

7.751  xlO’02 

^ laero 

1.681xlO+0° 

1.728  xlO+0° 

1.562  xlO+0° 

Table  5:  Computed  aerodynamic  forces,  moments,  and  power  for  prescribed  feath¬ 
ering  variation  using  the  4-term  kinematic  profile 

0%  Feathering  50%  Feathering  100%  Feathering 


CFx 

1.418xl0-°2 

5.116xl0-°3 

-3.927xl0-°3 

CFy 

-2.970xl0-°2 

-1.147xl0-°2 

4.607xl0-°3 

CFz 

1.352xlO-01 

1.506X10-01 

1.445  xlO"01 

CMx 

-8.513xl0-°2 

-8.921  xlO”02 

-8.582xl0-°2 

CMy 

-1.710xl0-°2 

-1.490  xlO”02 

-1.043xl0-°2 

CMz 

1.126xl0-°2 

9.512xl0-°3 

5.743  xlO”03 

CP 

8.360xl0-°2 

9.010xl0-°2 

9.640  xlO’02 

IJaero 

1.617xlO+0° 

1.672  xlO+0° 

1.499  xlO+0° 
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6. 5  Summary 

The  spanwise  feathering  variation  influences  the  orientation  of  the  wing  surface 
and  the  velocity  and  acceleration  profiles  of  the  wing  associated  with  the  feathering 
rotation.  This  spanwise  variation  in  the  kinematic  profile  influences  the  local  un¬ 
steady  aerodynamic  phenomena,  which  in  turn  influences  the  overall  aerodynamic 
performance. 

In  the  case  of  a  rigid  wing,  the  feathering  rotation  at  the  wing  root  induces 
a  non-circulatory  aerodynamic  phenomenon  that  produces  a  region  of  low  pressure 
on  the  “lower  wing  surface”  and  a  localized  set  of  vortical  flow  structures  near  the 
wing  root  during  pronation  and  supination.  At  the  same  point  in  the  wing  stroke, 
the  combination  of  circulatory  and  non-circulatory  aerodynamic  phenomena  at  the 
wing  tip  produce  a  region  of  high  pressure  on  the  “lower”  wing  surface  and  a  localized 
region  of  low  pressure  on  the  “upper”  wing  surface.  This  “reversed”  force  distribution 
decreases  the  net  aerodynamic  forces  produced  by  the  wing  during  the  rotational  phase 
of  motion.  The  localized  vortical  structures  at  the  wing  root  and  the  associated  region 
of  low  pressure  convect  below  the  wing  as  the  wing  changes  direction  of  motion.  This 
region  of  low  pressure  influences  the  pressure  distribution  along  the  wing  surface  as 
the  wing  sweeps  above  the  vortical  structures,  reducing  the  net  aerodynamic  forces 
produced  early  in  the  wing  stroke. 

In  addition  to  the  direct  influence  on  aerodynamic  performance,  the  local  pres¬ 
sure  distribution  and  vortical  structures  at  the  wing  root  of  the  rigid  wing  disrupt 
the  formation  of  the  leading  edge  vortex  on  the  inboard  part  of  the  wing  after  stroke 
reversal.  This  may  account  for  discrepancies  between  early  rigid  wing  simulations  and 
experimental  observations  of  the  flexible  natural  flier  regarding  the  continuity  of  the 
leading  edge  vortex  across  the  body  from  wing-tip  to  wing-tip  [21,92], 
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VII.  Analysis  of  Bending  Along  the  Elastic  Axis 

As  with  torsional  flexibility,  the  lack  of  a  good  characterization  of  the  Hawkmoth 
wing  structure  necessitated  the  use  of  assumed  bending  characteristics  rather  than 
the  use  of  a  structural  model  to  examine  the  influence  of  bending  on  the  unsteady 
aerodynamic  mechanisms.  Thus,  the  discussion  of  bending  influence  is  presented  as  a 
proof-of-concept  demonstration  and  used  to  identify  coupling  mechanisms  that  could 
influence  fluid-structure  interactions  of  interest  and  identify  future  research  directions. 

7.1  Prescribed  Bending  as  a  Function  of  Span  and  Time 

For  consistency  with  the  assumed  torsion  distribution,  the  first  bending  mode 
of  a  fixed-free  beam  was  used  to  define  the  spatial  bending  distribution.  The  first 
bending  mode  for  a  fixed-free  beam,  defined  by  Equation  (28)  for  the  coefficients 
a i  =  1.87510  and  b\  =  0.734096,  was  scaled  by  the  inertial  forces  due  to  the  prescribed 
wing  kinematics  to  approximate  the  dynamic  wing  response. 

r s(r)  =  cosh  ( arr/R )  —  cos  ( a\r/R )  —  b\  (sinh  ( apr/R )  —  sin  ( a\r/R ))  (28) 

The  time  history  of  bending  was  determined  using  the  inertial  forces  acting  on 
the  wing  during  the  wing  stroke,  which  depend  on  the  mass  distribution  of  the  wing 
and  the  kinematic  accelerations  in  sweep,  elevation  and  feathering.  By  co-locating 
the  elastic  axis  and  prescribed  feathering  axis  at  the  leading  edge  of  the  wing,  the 
influence  of  the  feathering  rotation  on  bending  may  be  neglected.  To  compute  the 
time-history  of  the  inertial  forces,  the  rotational  accelerations  due  to  feathering  and 
elevation  defined  in  the  inertial  reference  frame  were  transformed  into  the  wing  local 
frame.  This  transformation  provided  accelerations  normal  to  the  wing  surface  and  in 
the  chordwise  direction.  Lacking  a  well-defined  mass  or  stiffness  distribution  of  the 
Hawkmoth  wing,  the  time  history  of  the  accelerations  normal  to  the  wing  surface  and 
in  the  chordwise  direction  were  unit  normalized  and  scaled  to  achieve  a  reasonable 
maximum  tip  deflection.  Thus,  the  local  time-dependent  wing  deformations  were 
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defined  in  space  using  the  first  mode  of  a  uniform  fixed-free  beam  and  scaled  in  time 
by  the  rotational  accelerations  associated  with  the  kinematic  rotations. 

1.2  Results  of  the  Bending  with  Torsion  Simulations 

After  examining  the  available  footage  of  the  Hawkmoth  in  flight  and  considering 
the  smoothing  influences  of  the  torsionally  flexible  wings  on  the  unsteady  aerodynam¬ 
ics  at  the  wing  root,  the  bending  degree-of- freedom  was  superimposed  onto  the  50% 
torsion  simulation  described  in  Chapter  VI.  By  superimposing  the  bending  defor¬ 
mations  on  the  torsionally  flexible  wing,  illustrated  in  Figure  48,  the  bending  with 
torsion  simulations  provide  a  more  realistic  representation  of  the  flexible  Hawkmoth 
wing  than  would  be  achieved  by  superimposing  the  bending  deformations  onto  the 
rigid  wing  motion.  Two  peak  tip  deflections  were  considered  for  each  kinematic  profile 
representing  5%  and  10%  of  the  wing  length.  The  kinematic  accelerations  associated 
with  the  3-  and  4-term  kinematic  profiles  produce  two  time-dependent  tip  paths, 
shown  in  Figures  69  and  70. 


Figure  69:  Tip  deflection  history  for  the  3-term  kinematic  profile. 
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Figure  70:  Tip  deflection  history  for  the  4-term  kinematic  profile. 

Examination  of  the  force  history  of  the  bending  solutions  based  on  the  3-term 
kinematic  profile,  shown  in  Figure  71,  indicates  the  inclusion  of  bending  modes  in¬ 
creases  the  force  magnitude  at  the  end  of  each  half-stroke  while  decreasing  the  force 
magnitude  during  the  translational  phase  of  motion  in  the  middle  of  each  half  stroke. 
This  alters  the  3-term  force  history  from  a  single  aerodynamic  peak  during  each 
half-stroke  to  a  set  of  double  peaks.  As  would  be  expected,  the  increase  in  the  force 
magnitude  at  the  ends  of  the  wing  stroke  is  strongest  for  the  10%  tip  deflection.  How¬ 
ever,  due  to  the  wing  orientation  at  the  end  of  each  half-stroke,  the  increase  in  force 
magnitude  at  the  ends  of  the  half-stroke  has  a  greater  influence  on  the  ^-direction 
forces  than  the  z- direction  due  to  wing  orientation.  As  would  be  expected  for  first 
mode  beam  bending,  examination  of  the  spanwise  force  distribution  indicates  that 
the  largest  change  in  aerodynamic  forces  occcurs  near  the  wing  tip.  The  increase  in 
the  ^-direction  force  is  highlighted  in  Figure  72  at  t/r  =  8/8,  while  the  variation  in 
the  vertical  force  is  visible  in  Figure  73.  The  variation  in  the  tip  path  due  to  the 
bending  deformations  may  also  be  seen  Figure  73  The  influence  of  the  tip  deflections 
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Figure  71:  Aerodynamic  force  and  moment  histories  for  the  0%,  5%,  and  10% 

bending  simulations  using  the  3-term  kinematic  profile.  The  flier  is  oriented  in  a 
reference  system  with  the  x-axis  aligned  in  the  aft  direction,  the  z-axis  aligned  in  the 
vertical  direction,  and  the  y-axis  aligned  with  the  right  wing. 
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on  the  unsteady  aerodynamic  mechanisms  produced  by  the  3-term  kinematic  profile 
may  be  seen  as  variations  in  the  wake  structures  at  pronation  shown  in  Figure  74. 


Figure  72:  Comparison  of  the  spanwise  force  distribution  for  the  bending  simula¬ 

tions.  The  three  components  of  the  sectional  aerodynamic  force  are  shown  in  red  for 
the  ^-direction,  in  green  for  the  ^-direction,  and  in  blue  for  the  z- direction.  The  wing 
motion  is  based  on  the  3-term  kinematic  profile. 

The  prescribed  bending  deformations  based  on  the  4-term  kinematic  profile  in¬ 
fluence  the  force  history  in  a  similar  manner  as  the  3-term  profile,  increasing  the  force 
magnitude  at  the  ends  of  the  half-stroke  and  decreasing  the  force  magnitude  during 
the  translational  phase  of  motion,  as  shown  in  Figure  75.  Unlike  the  solutions  based 
on  the  3-term  kinematic  profile,  the  4-term  solutions  do  not  exhibit  the  same  dramatic 
increase  in  forces  at  the  beginning  the  half-stroke;  instead  the  deformations  increase 
the  amplitude  of  the  “shoulder”  discussed  for  the  4-term  feathering  solutions  in  the 
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Figure  73:  Comparison  of  the  spanwise  force  distribution  for  the  bending  simula¬ 

tions.  The  three  components  of  the  sectional  aerodynamic  force  are  shown  in  red  for 
the  ^-direction,  in  green  for  the  ^-direction,  and  in  blue  for  the  z- direct  ion.  The  wing 
motion  is  based  on  the  3-term  kinematic  profile. 
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Figure  74:  Iso-surface  of  vorticity  magnitude  for  the  0%  and  10%  Bending  simula¬ 
tions  based  on  the  3-term  kinematic  profiles  shown  during  pronation. 
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previous  chapter.  For  the  feathering  solutions  without  bending,  it  was  hypothesized 
that  the  shoulder  in  the  force  history  is  due  to  wing-wake  interactions,  while  the  early 
aerodynamic  peak  exhibited  by  the  4-term  kinematic  profile  was  due  to  the  circulatory 
aerodynamic  mechanisms  produced  when  the  peak  sweep  and  elevation  rates  occur 
in  phase  with  a  maximum  feathering  orientation,  ft  is  a  logical  extension  that  the 
bending  deformations  influence  the  wing-wake  interactions  by  increasing  the  effective 
tip  velocity  during  this  part  of  the  wing  stroke.  Wing-wake  interactions  could  account 
for  the  fluctuation  in  the  tip  forces  at  t/r  =  1/8  and  2/8  seen  in  Figures  76  and  77, 
as  well  as  the  variation  in  the  strength  and  coherency  of  the  tip  vortex  between  the 
0%  and  10%  bending  solutions  shown  in  Figure  80. 

The  response  late  in  the  half-stroke  of  the  4-term  kinematic  profile  is  dominated 
by  the  large  feathering  rotation  and  the  non- circulatory  aerodynamic  mechanisms  near 
the  wing  root.  The  dramatic  increase  in  force  magnitude  late  in  the  half-stroke  may 
be  due  to  non-circulatory  mechanisms  occurring  near  the  wing  tip  as  the  bending 
deformations  alter  the  effective  sweep  angle,  tip  velocity,  and  tip  acceleration.  This 
explanation  conld  account  for  the  variation  in  the  spanwise  force  distribution  seen 
in  Figures  78  and  78  at  t/r  =  3/8  and  t/r  =  4/8.  Examining  the  wake  structures 
produced  by  the  0%  and  10%  bending  solutions,  the  10%  bending  solution  exhibits  a 
significant  vortical  structure  shed  from  the  trailing  edge  of  the  wing  at  t/r  =  4/8  in 
Figure  80. 

Just  as  the  bending  deformations  influencing  the  wing-wake  interactions  by  in¬ 
creasing  the  tip  velocity  at  the  beginning  of  the  half-stroke,  the  deformations  exhibited 
by  the  4-term  kinematic  profile  reduce  the  translational  velocity  at  the  wing  tip  during 
the  middle  of  the  half-stroke  resulting  in  a  further  reduction  in  the  force  magnitude 
and  the  vertical  forces  shown  in  Figure  75.  The  variation  in  the  tip  path  due  to  the 
bending  deformations  may  be  seen  from  different  perspectives  in  spanwise  force  plots. 
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Figure  75:  Aerodynamic  force  and  moment  histories  for  the  0%,  5%,  and  10% 

bending  simulations  using  the  4-term  kinematic  profile.  The  flier  is  oriented  in  a 
reference  system  with  the  x-axis  aligned  in  the  aft  direction,  the  z-axis  aligned  in  the 
vertical  direction,  and  the  y-axis  aligned  with  the  right  wing. 
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Figure  76:  Comparison  of  the  spanwise  force  distribution  for  the  bending  simula¬ 

tions.  The  three  components  of  the  sectional  aerodynamic  force  are  shown  in  red  for 
the  ^-direction,  in  green  for  the  ^-direction,  and  in  blue  for  the  z- direct  ion.  The  wing 
motion  is  based  on  the  4-term  kinematic  profile. 
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t/r  =  1/8 


t/r  =  2/8 


Figure  77:  Comparison  of  the  spanwise  force  distribution  for  the  bending  simula¬ 

tions.  The  three  components  of  the  sectional  aerodynamic  force  are  shown  in  red  for 
the  ^-direction,  in  green  for  the  ^-direction,  and  in  blue  for  the  z- direction.  The  wing 
motion  is  based  on  the  4-term  kinematic  profile. 
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Figure  78:  Comparison  of  the  spanwise  force  distribution  for  the  bending  simula¬ 

tions.  The  three  components  of  the  sectional  aerodynamic  force  are  shown  in  red  for 
the  ^-direction,  in  green  for  the  ^-direction,  and  in  blue  for  the  z- direct  ion.  The  wing 
motion  is  based  on  the  4-term  kinematic  profile. 
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t/r  =  4/8 


Figure  79:  Comparison  of  the  spanwise  force  distribution  for  the  bending  simula¬ 

tions.  The  three  components  of  the  sectional  aerodynamic  force  are  shown  in  red  for 
the  ^-direction,  in  green  for  the  ^-direction,  and  in  blue  for  the  z- direction.  The  wing 
motion  is  based  on  the  4-term  kinematic  profile. 
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Figure  80:  Iso-surface  of  vorticity  magnitude  for  the  0%  and  10%  Bending  simula¬ 
tions  based  on  the  4-term  kinematic  profiles  shown  during  pronation. 
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7.3  Influence  on  Aerodynamic  Performance 

For  the  deformation  history  utilized  in  this  proof-of-concept  examination,  the 
mean  vertical  force  increases  with  increasing  tip  deflection,  shown  in  Table  6  for  the 
3-term  kinematic  profile  and  Table  7  for  the  4-term  kinematic  profile.  However,  this 
increase  in  lifting  capacity  was  tempered  by  an  increase  in  the  aerodynamic  power 
requirement.  For  both  of  the  kinematic  profiles,  the  torsionally  flexible  wing  with  no 
bending  represents  the  most  efficient  configuration. 


Table  6:  Computed  aerodynamic  forces,  moments,  and  power  for  prescribed  feath¬ 
ering  and  bending  using  the  3-term  kinematic  profile 


0% 

5% 

10% 

CFx 

3.638xl0_°3 

3.718x10-°° 

3.079x  10_°3 

CFy 

-1.625xlO-02 

-1.766xl0-°2 

-1.794xl0~°2 

CFz 

1.249  xlO”01 

1.285xlO_01 

1.376xlO"01 

CMx 

-7.608xl0-02 

-7.936  xlO”02 

-8.614xl0-°2 

CMy 

-1.280xl0~°2 

-1.350xl0-°2 

-1.314xl0-°2 

CMz 

4.001xl0~°3 

2.768  xlO”03 

7.654x  10~°4 

CP 

7.229x  10~°2 

7.717xl0-°2 

8.570xl0“°2 

Tjaero 

1.728x10+°° 

1.665x10+°° 

1.606x10+°° 

Table  7:  Computed  aerodynamic  forces,  moments,  and  power  for  prescribed  feath¬ 
ering  and  bending  using  the  4-term  kinematic  profile 


0% 

5% 

10% 

CFx 

5.117xl0-°3 

4.820  xl0"°3 

4.070xl0”°3 

CFy 

-1.147xl0-°2 

-1.417xl0-°2 

-1.814X10”02 

CFz 

1.507xl0-01 

1.552xlO_01 

1.619xlO_01 

CMx 

-8.922xl0~°2 

-9.280  xlO”02 

-9.811x  10“°2 

CMy 

-1.491X10"02 

-1.468xl0-°2 

-1.531x  10-02 

CMz 

9.512xl0~°3 

8.269  xlO”03 

6.694xl0-°3 

CP 

9.011xl0-°2 

9.692xlO-02 

1.060x  10~01 

Vaero 

1.672x10+°° 

1.601x10+°° 

1.527x10+°° 
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VIII.  Summary  and  Conclusions 

8.1  Summary  of  Research 

Navier-Stokes  based  simulation  was  utilized  to  understand  how  the  dynamic 
response  of  the  wing  structure  influences  the  aerodynamic  behavior  of  a  characteris¬ 
tic  flapping  wing  configuration.  The  use  of  Navier-Stokes  based  simulation  enabled 
detailed  examination  of  the  time-resolved  aerodynamic  forces,  unsteady  aerodynamic 
phenomena,  and  variations  in  the  time-resolved  wake  structures  for  the  large  ampli¬ 
tude  wing  rotations  utilized  by  natural  fliers.  The  influence  of  structural  deformation 
on  the  aerodynamic  response  was  examined  by  superimposing  torsional  and  bending 
modes  onto  the  rigid  wing  approximation  of  the  dynamic  wing  motions  utilized  by 
previous  numerical  investigations. 

8.2  Geometric  Configuration 

A  characteristic  configuration  based  on  the  Hawkmoth  Manduca  sexta  was  de¬ 
veloped  after  considering  the  influence  of  geometric  variations,  such  as  wing-body 
interaction,  wing  planform,  and  wing  length,  on  the  unsteady  aerodynamic  phenom¬ 
ena  and  aerodynamic  performance.  The  characteristic  configuration  developed  for  this 
study  utilized  a  single  wing  with  a  rectangular  planform  with  prescribed  wing  mo¬ 
tions  based  on  the  kinematic  parameters  documented  by  Willmott  and  Ellington  [6]. 
This  configuration  represented  a  minimum  geometric  configuration  that  reproduced 
the  characteristic  unsteady  aerodynamic  phenomena,  such  as  the  leading  edge  vortex 
and  wing-wake  interactions,  previously  documented  in  the  published  literature  for  the 
Hawkmoth  [15,16,19].  During  the  development  of  the  characteristic  configuration,  it 
was  observed  that: 

•  There  is  a  minimum  separation  distance  between  the  wing  and  body,  or  wing- 
body  gap,  that  is  required  in  order  to  resolve  the  wing  motion.  The  minimum 
wing-body  gap  depends  on  the  wing  planform,  the  origin  of  rotation,  and  the 
kinematic  parameters. 
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•  It  was  shown  that  the  aerodynamic  interactions  between  the  wing  and  body 
were  small  for  the  wing-body  gap  required  to  resolve  the  rigid  wing  motion  of 
the  Hawkmoth  planform  without  experiencing  contact  between  the  wing  and 
body  surfaces. 

•  Wing  planform  influences  the  amplitude  of  the  time-history  of  the  aerodynamic 
forces.  However,  the  character  of  the  force  history  and  wake  structures  were 
similar  for  the  Hawkmoth  and  rectangular  wing  planforms  considered.  This  in¬ 
dicates  that  the  Hawkmoth  and  rectangular  wing  planforms  both  induce  similar 
unsteady  aerodynamic  phenomena  for  the  same  dynamic  wing  motion. 

•  One  important  contribution  of  the  configuration  study  was  the  observation  that 
the  magnitude  of  the  aerodynamic  forces  produced  by  wings  of  varying  length, 
but  similar  planform  and  flapping  frequency,  scale  with  the  square  of  the  wing 
tip  velocity. 

8.3  Kinematic  Profile 

A  cursory  examination  of  the  Fourier  approximation  of  the  kinematic  profile 
highlighted  the  influence  of  several  kinematic  parameters  on  the  aerodynamic  response 
during  the  translational  and  feathering  phases  of  motion: 

•  Sweep  rate,  elevation  rate,  and  feathering  orientation  influence  the  aerodynamic 
forces  during  the  translational  phase  of  the  wing  motion  where  the  leading  edge 
vortex  and  the  circulatory  aerodynamic  mechanisms  play  an  important  role. 

•  The  feathering  acceleration  influences  the  non-circulatory  aerodynamic  phenom¬ 
ena  and  the  rate  of  change  in  the  bound  circulation  during  the  rotational  phase 
of  motion.  This  influences  the  transition  of  the  “upper”  and  “lower”  wing  sur¬ 
faces  as  the  suction  and  pressure  roles  alternate  at  the  end  of  each  half-stroke. 
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8-4  Influence  of  Feathering  Variation 

Local  structural  deformations  were  superimposed  onto  the  rigid  wing  motion  to 
examine  the  influence  of  structural  response  on  local  aerodynamic  response.  Lacking 
a  characterization  of  the  Hawkmoth  wing  structure,  documented  feathering  variations 
for  a  Hawkmoth  at  hover  were  approximated  using  the  first  torsion  mode  for  a  fixed- 
free  beam.  A  good  match  to  the  documented  feathering  variation  was  achieved  by 
reducing  the  amplitude  of  the  feathering  orientation  at  the  wing  root  while  maintain¬ 
ing  the  rigid  wing  feathering  orientation  near  the  wing  tip.  Wings  of  varying  torsional 
stiffness  were  approximated  by  varying  the  amplitude  of  the  feathering  orientation  at 
the  wing  root. 

•  The  feathering  rotation  dominates  the  motion  of  the  wing  near  the  wing  root. 

•  The  spanwise  variation  in  the  feathering  rotation  alters  the  velocity  and  accel¬ 
eration  profile  of  the  wing  away  from  the  elastic  axis. 

•  The  time-dependent  force  histories  and  the  mean  aerodynamic  force  and  power 
requirements  exhibit  a  strong  dependence  on  the  amplitude  of  the  span-wise 
feathering  variation. 

•  The  reduced  amplitude  of  the  feathering  acceleration  has  a  significant  influence 
on  the  local  non-circulatory  aerodynamic  phenomena  near  the  wing  root.  The 
change  in  the  non- circulatory  forces  has  a  direct  influence  on  the  local  pressure 
distribution  along  the  wing  length,  altering  the  net  aerodynamic  forces  during 
the  rotational  phase  of  motion. 

•  The  variation  in  the  non-circulatory  aerodynamic  phenomena  at  the  wing  root 
due  to  the  change  in  the  feathering  acceleration  influences  the  strength  of  the 
root  vortex.  This  change  in  the  root  vortex  alters  the  interaction  with  the 
circulatory  aerodynamic  phenomena  along  the  wing  length.  These  interactions 
in  turn  influence  the  circulatory  aerodynamic  phenomena  near  the  wing  tip  and 
the  non-circulatory  wing-wake  interactions,  and  have  a  direct  influence  on  the 
overall  aerodynamic  performance  at  hover. 
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•  The  feathering  variation  alters  the  wing  orientation  and  geometric  angle  of  at¬ 
tack  along  the  wing  span  during  the  translational  phase  of  motion.  This  change 
in  the  wing  orientation  influences  the  strength  of  the  leading  edge  vortex  and 
circulatory  aerodynamic  mechanisms  during  the  translational  phase  of  the  wing 
motion. 

For  the  kinematic  profiles  studied  in  detail,  the  solution  with  a  50%  reduction 
in  the  feathering  amplitude  at  the  wing  root  produced  the  greatest  mean  vertical 
force  and  represented  the  most  aero  dynamically  efficient  solution.  This  indicates  that 
wings  of  moderate  torsional  stiffness  may  outperform  both  highly  compliant  and  rigid 
wing  structures. 

8.5  Influence  of  Bending 

Local  deformations  based  on  the  first  bending  mode  of  a  fixed-free  beam  were 
superimposed  onto  the  torsionally  fiexible  wing  to  examine  the  influence  of  bending  on 
the  local  aerodynamic  phenomena.  The  time  history  of  the  bending  was  based  on  the 
inertial  forcing  of  the  wing  structure  due  to  the  kinematic  profile  and  the  amplitude 
of  the  bending  was  scaled  to  achieve  tip  displacements  of  5%  and  10%  of  the  wing 
length.  The  bending  deformations  influenced  the  aerodynamic  response  during  both 
the  translational  and  rotational  phases  of  motion. 

•  Incorporating  bending  deformations  increased  the  aerodynamic  force  magnitude 
during  the  rotational  phase  of  motion  while  decreasing  the  aerodynamic  force 
magnitude  during  the  translational  phase  of  motion. 

•  An  examination  of  the  vortical  wake  structures  suggests  that  the  force  variation 
early  in  the  half-stroke  are  due  to  enhanced  wing-wake  interactions. 

•  The  force  decrease  during  the  translational  phase  of  motion  corresponds  to  a 
reduction  in  tip  velocity. 

•  The  force  increase  at  the  end  of  the  half-stroke  is  due  to  additional  non-circulatory 
aerodynamic  influences  near  the  wing  tip. 
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•  The  influence  of  the  bending  deformations  on  overall  performance  was  to  in¬ 
crease  the  mean  vertical  force  while  decreasing  aerodynamic  efficiency. 

8. 6  Contributions 

Through  an  examination  of  the  influence  of  prescribed  structural  response  on  the 
unsteady  aerodynamic  phenomena,  this  work  provides  new  insight  into  the  influence 
of  structural  response  on  the  local  unsteady  aerodynamic  phenomena  and  overall 
aerodynamic  performance.  The  results  of  this  work  indicate  that  dynamic  structural 
response  plays  an  important  role  in  determining  the  aerodynamic  performance  of 
flapping  wing  fliers. 

In  the  case  of  a  torsionally  flexible  wing,  spanwise  variations  in  the  feathering 
orientation  influence  the  strength  and  timing  of  the  circulatory  and  non-circulatory 
unsteady  aerodynamic  phenomena  produced  by  the  dynamic  wing  motion.  These 
local  variations  in  the  unsteady  aerodynamic  phenomena  influence  the  evolution  of 
the  vortical  wake  structures,  surface  pressure  distributions,  and  the  directionality  and 
strength  of  the  net  aerodynamic  forces. 

In  the  case  of  the  flexible  wing  with  bending  degrees  of  freedom,  the  circula¬ 
tory  and  non-circulatory  unsteady  aerodynamic  phenomena  were  influenced  by  the 
variation  in  the  span-wise  velocity  profile,  particularly  during  the  rotational  phase  of 
motion.  The  sensitivity  of  the  aerodynamic  performance  to  bending  degrees  of  free¬ 
dom  was  observed  for  tip  deflections  of  5%  and  10%  of  the  wing  length.  This  indicates 
that  spanwise  bending  may  still  play  an  important  role  in  determining  aerodynamic 
performance  for  wings  with  pronounced  variations  of  the  effective  bending  stiffness  in 
the  spanwise  and  chordwise  directions. 

This  work  also  highlights  the  sensitivity  of  the  global  wake  structure  to  small 
variations  in  the  local  aerodynamic  phenomena.  This  sensitivity  was  observed  through 
qualitative  comparisons  of  the  vortical  wake  structures  and  quantitative  comparisons 
of  the  net  aerodynamic  force  histories.  This  sensitivity  of  the  aerodynamic  forces  and 
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wake  structure  to  small  variations  in  the  unsteady  aerodynamic  phenomena  high¬ 
lights  the  need  to  adequately  resolve  the  underlying  unsteady  aerodynamic  phenom¬ 
ena  within  future  design  tools  that  consider  the  influence  of  structural  response  and 
kinematic  variations  on  aerodynamic  performance. 

8.7  Recommendations  for  Future  Research 

The  present  work  provides  a  benchmark  study  examining  the  influence  of  struc¬ 
tural  response  on  aerodynamic  performance  using  Navier-Stokes  based  simulation.  By 
extending  the  simulation  tools  developed  for  this  study,  the  role  of  multi-physics  cou¬ 
pling  in  determining  system  performance  may  be  explored  by  examining  the  following 
questions: 

•  Do  aerodynamic  forces  or  inertial  clastic  forces  dominate  the  structural  re¬ 
sponse? 

•  How  does  the  phase  relationship  between  the  aerodynamic  and  inertial-elastic 
forces  influence  the  structural  response? 

•  Is  the  feathering  rotation  due  to  passive  structural  response  or  active  actuation 
of  the  wing  structure? 

•  Are  higher-order  modes  present  in  the  structural  response  of  the  coupled  system? 

—  If  so,  what  influence  do  higher  order  modes  have  on  the  unsteady  aerody¬ 
namic  phenomena? 

•  How  does  fluid-structure-kinematic  coupling  influence  agile  maneuvers? 

—  What  is  important  to  the  transient  response  of  the  system? 

—  What  is  the  phase  lag  between  change  in  wing  actuation  and  change  in  the 
force  history? 

—  If  there  is  a  phase  lag,  is  it  due  to  structural  response  or  aerodynamic 
response? 
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In  addition  to  the  role  of  multi-physics  coupling  on  aerodynamic  performance, 
the  sensitivity  of  the  unsteady  aerodynamic  phenomena  to  prescribed  changes  in  the 
geometric  configuration,  kinematic  profile,  and  structural  response,  may  be  explored 
using  the  current  set  of  simulation  tools: 

•  What  is  the  influence  of  Reynolds  number  on  the  unsteady  aerodynamic  phe¬ 
nomena,  and  how  should  Reynolds  number  be  defined  for  hovering  configura¬ 
tions? 

•  What  is  the  role  of  turbulence  in  insect  flight,  and  how  important  is  laminar- 
turbulent  transition  for  capturing  the  near-body  and  off-body  wake  structures? 

•  How  does  the  choice  of  wing  planform  influence  aerodynamic  performance? 

—  Is  there  an  ideal  spanwise  force  distribution  for  efficiency  and  power? 

—  What  influence  does  area  distribution  and  taper  ratio  have  on  the  aerody¬ 
namic  phenomena? 

—  How  does  tip  shape  influence  the  interaction  between  the  leading  edge  and 
tip  vortices? 

•  What  is  the  influence  of  the  origin  of  rotation  on  aerodynamic  performance? 

—  How  does  the  location  of  the  feathering  axis  influence  the  leading  edge 
vortex? 

—  How  does  the  location  of  the  feathering  axis  influence  the  non-circulatory 
phenomena  during  pronation  and  supination? 
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